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ABSTRACT 


Application of classic triangulation methods will allow 
the location of a radar to be determined by passive sensors. 
Through the use of modern digital signal processing techniques 
this estimate can be made in a simpler fashion using a 
conventional receiver. 

In this thesis a technique is developed for time 
difference of arrival (TDOA) estimation using a frequency 
domain based correlation detector driven by an envelope 
detector. Time lag boundaries are defined on the output of the 
correlator. A fixed detection threshold is calculated to 
permit constant false alarm rate (CFAR) detection. The 
performance of the correlation detector is plotted as a 
receiver operating characteristic (ROC) curve as a function of 
signal to noise ratio (SNR). An interactive MATLAB software 
program is provided to perform either spectral domain or time 
domain based correlation. 

Spectral domain based correlation uses the Fast Fourier 
Transform (FFT). Implicit with the use of the FFT are finite 
arithmetic internal processing errors which are modeled as 
independent uncorrelated noise sources. A method is presented 


to account for SNR degradation at the output of the FFT. 
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I. INTRODUCTION 


A. BACKGROUND 

Interception of radar transmissions by passive sensors is 
one of the responsibilities of electronic intelligence. 
Analysis of the signal of a radar can provide an estimate of 
the distance and direction of an emitter relative to an 
intercepting receiver. Knowing the location is an important 
consideration in the evaluation of strategic and tactical 
capabilities of the radar. 

Triangulation is a traditional technique where the angle 
of arrival of a radar signal at two spatially separated 
receivers provides an estimate of the transmitters location. 
With the advent of the computer, digital signal processing 
algorithms can be implemented that provide emitter range 
information and hence allow localization in a sequential 
sense. One of these processing algorithms is called the time 


difference of arrival (TDOA) method. 


B. TIME DIFFERENCE OF ARRIVAL 

The TDOA algorithm relies on two bistatic passive sensors 
attempting to receive the transmission of an active radar. 
Reception by the two sensors is preferentially done in the 


main lobe of the transmitting radar, but due to the 


nonsynchronous nature of the task, reception through one of 
the sidelobes of the transmitting radar is expected with a 
resulting decrease in received signal power. 

Given a transmitted radar pulse f(t), two physically 
separated passive sensors will most likely receive this pulse 
at times t, and t,tt. The time difference t is proportional to 
their differential distance relative to the emitter. This time 
difference of arrival can be used in the estimation of the 
differential range of the emitter to the receivers. The 
reception of the first pulse in a pulse burst by the sensors 
is paramount in performing the TDOA estimate. 

The concept of measuring the TDOA of a radar pulse between 
two separated receivers can be modeled by a hyperbola. The two 
sensors whose positions are known are located at the two focal 
points of the hyperbola. A transmitting radar is located in 
the area surrounding the two sensors. Figure 1 shows this 
model. 

A hyperbola is the locus of points in which the distance 
from any one point on the hyperbola to one focus differs by a 
constant amount from the distance to the other focus. This 
differential distance is proportional to a difference in 
arrival time t, of a radar pulse as detected between the two 
sensors. Once T is measured the hyperbola can be drawn. If the 
sensors (moving aircraft or nongeostationary satellites) 
repeat the differential arrival time measurements at several 
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time intervals, a series of hyperbolas will be generated, 
fixing the location of the transmitting radar.  Alternately, 
the angle of arrival at either sensor can be used to define 
the position of the transmitter on the hyperbola. A good 
introduction into the use of the hyperbola in navigation is 


found in [Ref. 1:p. 27). 


Range (meters) 


Transmitting Radar 


Bl 


Sensor 1 Sensor 2 Range (meters) 


(focal point 1) (foca! point 2) 


Hyperbola for S 1 ms. 
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Figure l. Time difference of arrival model. 


This discussion assumes that the radar signal received by 
the two sensors is induced by the direct wave from the 
transmitter. The effects of refraction and reflection on 


propagating radar waves are not covered here. 


The purpose of this thesis is to develop and test an 
algorithm to estimate the differential arrival time t, of a 
pulsed radar signal collected by two passive sensors. If the 
pulse burst is collected in the time domain, time correlation 
can be used to generate a TDOA estimate. The Fast Fourier 
Transform (FFT) can also be used to perform the correlation of 
the two received signals (i.e., fast correlation). 

Additive noise n(t), superimposed onto a transmitted radar 
signal f(t) by the environment and the radar receiver, creates 


a composite signal s(t) at the receiver 


st) = Eier (1) 


The use of the Fourier Transform for TDOA estimation is 
attractive because of the processing gain (PG) a radar signal 
receives in being transformed from the time domain to the 
frequency domain during the detection process. FFT processing 


gain for real valued signals can be approximated by 


PG (dB) = (log, [number points in data record] -1)'3. (2) 


See for example [Ref. 2:p. 33]. The signal to noise ratio 
(SNR) at the output of a FFT (SNR,,) is a function of the PG 
and the SNR at the input to the FFT (SNR). This relation is 


given by 


SNRo,, " SNR,, + PG . (3) 


The corruptive effects of additive noise are reduced by the 
transformation, allowing the frequency domain detection of the 


Signal. 


C. RADAR TYPES 

The TDOA algorithm is primarily applied to pulsed radars. 
The continuous wave (CW) radar provides target radial velocity 
through the Doppler shift of its return. The CW radar does not 
use a pulsed transmission. Therefore, the TDOA algorithm 
cannot not be directly applied to the reception of CW signals. 

Pulsed radars can be divided into two broad categories, 
ordinary low pulse repetition frequency (PRF) pulsed radars 
and the pulse Doppler radars. 

Low PREF radars can yield unambiguous range information, 
while their radar returns do not give target velocity 
information directly. Generally because of the long duty cycle 
of the low PRF radar, the probability of receiving all radar 
pulses in a pulse burst by two passive sensors is great. The 
TDOA algorithm can then estimate the differential distance to 
the emitter using the low PRF burst. 

The pulse Doppler radar uses coherent transmission and 
reception with a moderately high PRF [Ref. 3:p. 17.1]. It 
gives both range and radial velocity information about a 
target, though not with the same accuracy of the ordinary 
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pulsed radar or CW radar respectively. In the pulse Doppler 
radar, as in all non-bistatic radars, the receiver must be 
turned off while the transmitter is on. Because of this, a 
range ambiguity exists. Echoes having delay times equal to an 
integral multiple of the PRF will be undetected. Target echoes 
that remain within this blanked time interval will remain 
undetected. A second source of range ambiguity exists. For a 
target located at a distance r, from the tracking radar, all 


other targets in the main beam of the radar a distance 


I, T4; IL, A E 
(4) 
where k - integer , 
will appear to be approximately the same distance from the 
radar. 

Pulse Doppler radars can be divided into a low, medium and 
high PRF class. Both the low and medium PRF pulse Doppler 
radar transmissions can be passively intercepted and a 
relatively simple TDOA estimate can be formed. 

Reception of each pulse in a pulse burst from a high PRF 
pulse Doppler radar by two independent passive sensors can be 
difficult. If either sensor misses the first pulse, the TDOA 
estimate will be in error. This is true for any radar, but has 
a higher likelihood of occurring in the high PRF radars due to 
the shorter pulse widths used. Should the pulse burst be coded 
so that a nonuniform time pulse sequence is generated (i.e., 
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staggered PRF), then missed pulses by either sensor will not 


seriously degrade the validity of the TDOA estimate. 


D. FAST FOORIER TRANSFORM OUTPUT SIGNAL TO NOISE RATIO 

The Fourier transform is a computationally intensive 
mathematical operation. When implemented on an FFT processor 
using fixed point arithmetic, internal processing errors are 
generated. These errors are a function of the transform size, 
and the number of bits used internally in the FFT processor. 

The primary building block of the FFT is the butterfly 
algorithm. There are three dominant processing errors that 
occur within the butterfly. They are scaling, truncation and 
trigonometric errors. These can be modeled as independent, 
uncorrelated noise sources. The cumulative effect of these 
noise sources is to reduce the potential SNR at the output of 
the butterfly, and therefore, the FFT processor. 

This thesis will also examine the effects of noise caused 
by FFT processing, and will present a method to account for 
the degradation of the output SNR due to finite arithmetic. 
Certain FFT implementations will keep the signal energy larger 


with respect to the three noise power sources than others. 


II. TIME DOMAIN BASED CORRELATION 


A. BACKGROUND 
The correlation function R(T) measures the degree of 
similarity between two signals s(t) and g(t) and is described 


by 


ce 


AS: jc 
limit,.. ——| s(t)g(s + t)dt (5) 


where g(t + t) g(t) displaced by the shift x 


Correlation can be performed in a radar receiver between 
a target return corrupted by additive noise, and a replica of 
the transmitted signal maintained by the radar. This replica 
is designed into the frequency response of the matched filter 
of the radar receiver. [Ref. 4:p. 373]. 

Correlating two signals produces an output that does not 
resemble either of the two input signals. The shape of the 
received waveform is destroyed during the correlation process. 
The useful item of correlation is an amplitude peak, where the 
location along the time axis of this peak indicates the amount 
of time that one signal lags the other. 

If the correlation process occurs between a signal and a 
delayed replica of itself, the process is called 
autocorrelation. If two dissimilar signals are correlated, the 
process is called crosscorrelation. 
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B.  AUTOCORRELATION 

Autocorrelation is the process by which two identical 
signals, one a delayed replica of the other, are correlated. 
The autocorrelation function of a discrete time sequence x(i) 


is defined by [Ref. 5:p. 556] 


N-1-|1| 


Eee ex (i) ase + 1) 
1-0 (6) 
where 1 = shift operator 
N = number of data points in x(n) 


This function is not scaled with respect to the shift operator 
nor the number of data points. Figure 2(a) shows the effect of 
autocorrelation on a burst of 13 pulses of unit amplitude. 
Normalized autocorrelation produces a triangular envelope 
waveform of height one at a time lag of zero. The location of 
this peak indicates that the time shift for maximum overlap 
between the pulse burst and a shifted replica of itself is 
zero. 

For illustration, a Barker coded sequence of pulses of 
length 13 is also autocorrelated in Figure 2(b). The Barker 
code is given by +++++--++-+-+, where + and - could represent 
O and KT radians carrier phase shift, respectively. Barker 
coding of a pulse sequence constructively modifies the 
correlation output to magnify the time lag peak. The zero lag 
point has the greatest amplitude.  Adjacent peaks are 


significantly reduced in amplitude. Because of this reduction, 


uniform autocorrelation, 50% duty cycle : Barker code autocorrelation, 50% duty cycle. 


amplitude 
amplitude 


Hu | 
Al Hl | 





Figure 2. (a) Normalized autocorrelation of 13 equal 
amplitude pulses. (b) Normalized autocorrelation of 13 Barker 
coded pulses. 


if noise is added to the received signal, a Barker coded radar 
pulse burst will have a correlation peak that is not as easily 
confused with adjacent peaks. 

Now consider the autocorrelation of a random process. An 
independent,  identically distributed zero mean Gaussian 
sequence with variance of one is shown in Figure 3. The 
autocorrelation of this sequence produces a peak at zero time 
lag. There is no strong correlation of the Gaussian noise 
except at zero time lag. The autocorrelation of Rayleigh 


distributed noise is shown in Figure 4. 
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Gaussian distributed noise normalized Gaussian autocorrelation 


9 
= 
E 
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E 
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amplitude 





Figure 3. (a) Gaussian distributed noise sequence. 
(b) Normalized autocorrelation of Gaussian noise. 


Rayleigh distributed noise normalized Rayleigh autocorrelation 


amplitude 





Figure 4. (a) Rayleigh distributed noise sequence. 
(b) Normalized autocorrelation of Rayleigh noise. 
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Rayleigh distributed noise has a mean (dc) value of 


ol = 
2 (7) 


standard deviation 


mean value 


where o 
The autocorrelation function of Rayleigh noise has two 
components. A triangular envelope caused by the dc component 
of the Rayleigh noise, and a weighted delta function at the 
apex of the triangle which indicates the time shift of maximum 
overlap of the two sequences. Clearly maximum overlap occurs 


at a time shift of zero for any autocorrelation function. 


C.  CROSSCORRELATION 
Crosscorrelation is the process wherein two different 
signals are correlated. The crosscorrelation function of a 


discrete time sequence is defined by 


N-1-|1| 


is 2. x(i)y(i + 1) (8) 


where x(i) and y(i) are two different sequences. 
This function is not scaled with respect to the shift operator 
nor the number of data points. Figure 5 shows the 
crosscorrelation of two uncorrelated Gaussian noise sequences 
and two uncorrelated Rayleigh noise sequences. Clearly, no 
dominant amplitude peak is visible in either crosscorrelation 
plot. Lack of a dominant peak indicates there is no strong 


correlation between the two sets of noise sequences. 


E- 


normalized Gaussian crosscorrelation normalized Rayleigh crosscorrelation 





Figure 5. (a) Normalized Gaussian noise crosscorrelation. 
(b) Normalized Rayleigh noise crosscorrelation. 


D. CORRELATION COEFFICIENT 

A measure of the linear dependence between any two zero 
mean stationary time functions x(t) and y(t) is the 
correlation coefficient function and is defined [Ref. 6:p. 48] 


as 


Pxy (1) = — Iom _ 


(EO ON (9) 


where |p,, (y) | s 1 for allt. 


This function is the crosscorrelation function normalized by 
the square root of the maximum values of the individual 
autocorrelation functions. The normalizing factor is not a lag 


dependent quantity. 
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Given two time series x(t) and y(t) where a is any 


positive number 


Py [TY = 1 if XC E 
Pa (T) = -1 if x( OT- Ta (10) 
p, U) * O if x(t) and y(t) ere neem we 2 


Another way of defining a normalized correlation 


coefficient for a discrete time series of finite duration is 


given by 
N-1-|2| 
y x) y Em 
Py DY. A mů (11) 
N-1 N-1-|1| 


Soxa S 
0 1=0 


E. SQUARE-LAW DETECTION AND CORRELATION 
If x(t) is a real valued function of time, the Hilbert 


transform of x(t) is defined by [Ref. 6:p. 484] as 


X(C) 


HE) 


=f x(u) zi 
-.Tt(t€ - u) (12) 


ak 
x(t) * Smet 


where * = convolution operation 


The output of a square-law envelope detector u(t), can be 


described by 
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E qe t) 


(13) 


real valued time function 
input to square-law detector 


where x(t) 


The correlation coefficient function of two time functions 
x(t) and y(t) is defined in Equation 9. The Hilbert transform 


of the correlation function is defined as 


B, (1) =H [p,,(t)] 


= Ry (T) 


Qu hy, (14) 


standard deviation of x and y 
respectively 


where o, and 0, 


The correlation coefficient for two square-law envelope 


detected signals, p,,(t), is defined by [Ref. 6:p. 512] as 


pue p.c: D. (v) 
(15) 
where u(t) = square-law detected signal of x(t) 
v(t) = square-law detected signal of y(t) 


The quantity p,,(T) produces a sharper correlation peak than 
Pz (T). This sharpening of the correlation function peak output 
can aid in locating the TDOA correlation peak. An example of 


Pott P(T) and Pu (T) is plotted in Figure 6. 


F. ENVELOPE CORRELATION AND RELATED STATISTICS 

In this thesis, time domain based correlation will be 
performed on the output of an envelope detector. An envelope 
detector can be easily implemented in hardware using a diode. 
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normalized correlation coefficient normalized correlation coefficient and Hilbert 
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Figure 6. (a) Normalized correlation coefficient. (b) Hilbert 
transform of the normalized correlation coefficient (dotted 
line}. 
To derive the expected value and variance of the 
crosscorrelation function, two cases must be considered. These 
cases are noise only present, and signal plus noise present in 
both channels of the correlator. Under the noise only case, 
the statistics for the correlator output are easily derived. 
We consider two real, independent identically distributed, 
zero mean noise series at the input to the correlator. Each 
series has zero mean and variance 6! . It is shown in Appendix 
A, that the expected value of the output of the 


crosscorrelation functions z DE 


(1)] =0 . (16) 
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The variance of the output of the crosscorrelation function 


m (S 


2 = 4 
oc ES (17) 


where N - number of data points 


Since equal variances are assumed for both time series, 0° is 


the square of the variance of either noise series. 


G. BLOCK CORRELATION 

Two time sequences can be correlated by any of three 
different methods. These will be discussed below. 

The first method blocks the length of each input sequence. 
Correlation of these blocks will provide an output with a 
variance given by Equation 17. The linear dependence of the 
correlation variance on (N-|@|) is shown in Figure 5 for zero 
mean Gaussian noise. Figure 5 shows the correlated output at 
(=0 (zero time lag) is scaled by (N-|0])=N. Correlation output 
at a time lag of +N are scaled by (N-|N|)=0. The linear 
dependance of the correlation output variance on the time lag 
does not allow a constant detection threshold, which is needed 
to automatically determine the location of the correlation 
peak. 

A second method of correlation blocks a fixed length 
sequence x(n) and correlates it with an arbitrarily long 


sequence y(n). The correlation output for this algorithm has 
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a variance of No‘. Clearly, no scaling of the correjasısn 
output variance as a function of time lag would occur. 

A third method is scaling applied to the output of the 
first technique. The correlator output is multiplied by a 
weighting function, 1/(N-|@|)?/? to correct for the lag 
dependency of the variance. The selection of a constant 


detection threshold will be addressed in Chapter V. 
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III. FREQUENCY DOMAIN BASED CORRELATION 


A. FOURIER TRANSFORM 


Frequency domain techniques can be used to perform fast 
correlation rather than the time domain techniques previously 
discussed. The Fourier transform translates the discrete time 
sequence x(n) into the frequency domain 

For a finite duration time series x(n), the discrete 


Fourier transform (DFT) is defined by the pair of equations 


[Ref 5:p. 100] 


N-1 en 
-J = EK 
X(k) = x(n) e NU" 0 « k « N-1 
n-0 
(18) 
N-1 j( 2) 
x(n) = 2 x(k) e! "".osnsmN-i. 
N k= 


Equation 18 is used to transform the N discrete sample values 
into N frequency domain coefficients. 

The following discussion is based on [Ref. 5:p. 286]. The 
DFT is a computationally intensive algorithm if implemented 
directly. The direct evaluation of the DFT requires N’ complex 
multiplications. The time required for the evaluation of a DFT 
is proportional to N. 

Algorithms have been designed that take advantage of the 
symmetry and periodicity of the DFT, reducing the amount of 
computation required to solve the DFT. Any of these algorithms 


are referred to as the FFT. 
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Generally, the FFT requires (N/2) 1og;N complex 
multiplications. This smaller computational workload becomes 
significant as N, the number of data points in the time 


series, becomes larger. 


B. NORMALIZED FREQUENCY DOMAIN CORRELATION 

In some radar receivers the processing of pulse sequences 
occurs in the frequency domain. Therefore, a frequency 
representation of the time domain pulses may already exist. 
The frequency coefficients do not necessarily need to be 
transformed back into the time domain to perform the 
crosscorrelation. 

Frequency domain data can be  crosscorrelated by 
multiplying one set of coefficients by the complex conjugate 
of the second set of coefficients. The product is translated 
back to the time domain through the use of the inverse FFT. 
Frequency domain based normalized  crosscorrelation is 


described by 
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0 otherwise 


where X(k) = FFT of the series x(n) 
Y(k) = FFT of the series y(n) 
* = conjugation 
po tousenormsbdcedusemnenatdonetunetdon 


TO avoid circular correlation, both time domain pulse 


sequences are zero padded to N, data points 


T EE 


where N, = number of data points in x(n) (20) 
N- number of data points in y(n). 


To take advantage of the processing speed of FFT algorithms, 
N, is made a power of two by increasing its length with 


additional zero padding 


N = 27 ZN ee 
Z X y (21) 


where m is an integer. 


Prior to any zero padding, the dc component is removed to 
Create zero mean pulse sequences. Both sequences, each of 
length N, are then transformed to the frequency domain through 
the FFT and processed using the frequency domain correlation 
technique. A MATLAB implementation of frequency domain 


crosscorrelation is given in Appendix B (FreqCorr7.m). 


ZI 


C. TIME DIFFERENCE OF ARRIVAL (FOURIER DOMAIN) 

The normalized crosscorrelation function implemented in 
the frequency domain (Equation 19) is used to estimate the 
TDOA. The crosscorrelation of two similar signals will produce 
a dominant peak in the output of the correlator. This peak is 
used as an estimate of the TDOA between the two signals. For 
two signals that have not been corrupted by noise, the 
location in time of this peak is easily determined. Pulsed 
signals that have been distorted by noise produce a noisy 
correlated output. The output noise can mask the TDOA peak in 
low SNR conditions. Consequently, the probability of selecting 
a noise sample instead of the true TDOA peak increases with 
decreasing SNR. 

Figures 7,8 and 9 show the frequency domain based 
crosscorrelation of two pulse sequences as a function of 
decreasing SNR. A reception is simulated with a pulse burst 
demodulated by one receiver fed into one channel, and a pulse 
burst demodulated by a second receiver fed into the second 
channel of the correlator. Pulse burst one has a transmission 
delay of zero. Pulse burst two has a delay of 50 time units. 
Therefore, the total TDOA is 50 time units. The TDOA peak is 
clearly seen at the high SNR of +7dB and +3dB. At OdB SNR, 
noise peaks exceed in amplitude the true TDOA peak causing an 


incorrect estimate if one was selected. 
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Figure 7. Normalized frequency domain based crosscorrelation. 
SNR +7dB. 
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Figure 8. Normalized frequency domain based crosscorrelation. 
SNR +3dB. 


Visually estimating the location of the peak with the 
maximum amplitude in the output of the correlator is a easy 
task. The decision rule states that the point with the maximum 
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Figure 9. Normalized frequency domain based crosscorrelation. 
SNR OdB. 


amplitude is selected as the TDOA estimate, provided the 
amplitude exceeds a predetermined threshold. 

When many blocks of data must be processed, visual 
inspection of the correlation output is not possible. Chapter 
V discusses the design of a threshold for automated data 


processing. 
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IV. NOISE AND SIGNALS IN RECEIVER 


A. NOISE AND SIGNALS IN THE RADAR DETECTOR 

The purpose of a radar receiver is to process and detect 
the reflected energy from a target being tracked by the radar. 
An I/Q demodulator is assumed in the radar receiver, located 
after the last intermediate frequency (IF) amplifier and prior 
to the video signal display. Data to be processed for TDOA can 
be obtained from three locations on the I/Q demodulator. These 
locations include the envelope output, the envelope squared 
output and the I/Q channel data. 

White Gaussian distributed noise is assumed to be added to 
radar pulses during their transmission through the atmosphere 
and reception by the receiver. This noise contributes to the 
degradation of the SNR and is caused by fluctuations in the 
target radar cross-section and the loss of signal energy due 
to the propagation distance between the target and the 
receiver. 

At typical radar frequencies the front end amplifier of 
the radar receiver contributes to this noise power. The noise 


power for a thermally noise limited receiver is defined as N,, 
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No = Ki Gee 
where k = Boltzman’s constant 
B = IF bandwidth 
(22) 
T ES yr AA 


Tą = receive antenna temperature 
T, - effective noise temperature 
of the receiver 


The total noise is modeled as zero mean Gaussian distributed 


noise with a probability density function (pdf) of 


-x^ 


dl (23) 


poc a 5 


y2 "0? 
where O^ is the noise variance due to all noise sources (i.e., 
receiver front end, transmission noise). 

For an operational radar, either a pulsed signal plus 
Gaussian noise or Gaussian distributed noise alone is assumed 
at the input of the I/Q demodulator. The probability 
distribution of the demodulated signal is dependant on which 
of the three outputs of the I/Q demodulator is used. 

1.  Coherent detection 

If the exact frequency and phase of the received pulse 
burst is known, and matches the frequency and phase of the 
local oscillator in the receiver, then the signal can be 
coherently detected.  Coherent detection simplifies the 


probability description of the signal and noise at the output 
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of the I and Q channels. The I and Q sequences have an 


independent Gaussian distribution, and can be described as 


T Que NC Um (nieto ene, en (24) 
where m;(n), mg(n) * I and Q mean values respectively 
o? = noise variance due to all noise sources. 


Both the I or Q channel can be processed using the correlation 
algorithm followed by threshold detection. 

Due to the passive nature of signal reception and 
demodulation, the exact frequency and phase of the received 
pulse burst are not available at the receiver. Therefore, 
coherent detection is not feasible for TDOA estimation. 

2. Envelope Squared Detection 

For a zero mean Gaussian distributed noise input to 
the I/Q demodulator, I? and Q? each have a chi-squared 
distribution with one degree of freedom. Their pdf's can be 


described by 


By) 2o for y z 0 


= 0 for y TO (25) 


where 0° = variance of I and Q 
y= I oro. 


The envelope squared output is the sum of I? and Q? and is 
therefore chi-squared with two degrees of freedom. This pdf 


can be described by 
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E e 2/20? 
p,(z) = — TORE 
20 (26) 


where z = I? + Q? 


Envelope squared detection is attractive because the Hilbert 
transform can be used to enhance the TDOA peak in the 
correlator's output. Unfortunately in a low SNR environment 
the noise power will also increase as the noise contribution 
is squared. 
3. Envelope Detection 

A common form of receiver demodulation at IF 
frequencies is the envelope detector. For all simulations in 
this thesis, envelope detection will be assumed. The output of 
the envelope detector is the modulation envelope of the 
received pulse burst. For Gaussian distributed noise at the 
input to the I/Q demodulator, the noise at the envelope output 


is Rayleigh distributed with pdf given by 


u? 


=e 20° Uaz 0 
g? 
0, otherwise (27) 


To ur) 


where gł = variance of Gaussianmelse 


Vz = VI? + o? 


= 
if 


For a pulsed sinusoid plus Gaussian noise at the input to the 
I/Q demodulator, the pdf of the envelope at the output is 


Rician distributed. The Rician pdf is given by 
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where A = amplitude of received sinusoid 
0^ - variance of Gaussian noise 


2 
where I,(Z) à a [ " gzcost8) ap 
2T Jo 


(29) 
and I,(z) = modified Bessel function of the 
first kind of zero order 


Figure 10 shows a block diagram of the I/Q demodulator and the 


various pdfs in the radar receiver. 


B. ENVELOPE CORRELATION OF RAYLEIGH NOISE 

Rayleigh distributed noise is produced at the output of 
the envelope detector under noise only conditions. This noise 
is then processed by the correlator during the TDOA estimate. 
In all follow on discussions and simulations, the envelope 
processing scheme is used. 

The autocorrelation of Rayleigh noise as discussed earlier 
produces a large triangular bias in the output which decreases 
as the time lag increases. This bias is produced by the dc 
component of the Rayleigh distributed noise. 

This triangular function contains no information relative 
to the TDOA estimate. The large numerical values in the bias 
can potentially limit the dynamic range in the FFT based 
correlation. A constant, detection threshold cannot be 
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Figure 10. I/Q demodulator in the radar receiver. 


implemented at the output of the correlator with this 
triangular function present. For these three reasons, it is 
advantageous to force the Rayleigh or Rician distribution to 
have a zero mean. This will remove the triangular bias formed 
by the correlation process. This modification of the received 
pulse burst does not degrade the TDOA estimate. For all 
realizations (i.e., blocks of data), the dc component of the 


Rayleigh noise will be removed. 
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1. Shifted Rayleigh Noise 
Forming a Rayleigh distributed time series x(n) and 
subtracting its expected value produces a shifted Rayleigh 


sequence. This algorithm is described by 


x(n) =x(n)-o0 x 
: (30) 
where x,(n) » shifted Rayleigh noise sequence 
c = Gaussian noise standard deviation 


The shifted Rayleigh time series has a zero mean. The new time 
series and the corresponding time domain based autocorrelation 


function are shown in Figure 11. 


shifted Rayleigh time series normalized shifted Rayleigh autocorrelation 





Figure il. (a) Shifted Rayleigh time series. (b) Normalized 
shifted Rayleigh time domain based autocorrelation. 


Histograms of a Rayleigh pdf and a shifted Rayleigh 
pdf are shown in Figure 12. Clearly, the effect of removing 


the dc term from the Rayleigh distributed time series creates 
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Rayleigh pdf shifted Rayleigh pdf 





Figure 12. (a) Rayleigh pdf. (b) Shifted Rayleigh pdf. 


an autocorrelated output that does not have the limiting 
triangular bias. The impulse at zero time lag dominates the 
correlation output. A constant detection threshold, at least 
over a small number of range bins, could be implemented with 
this type of correlation output. 


2. Frequency Domain Based Correlation of Shifted Rayleigh 
Noise 


The frequency domain based autocorrelation of x,(n) is 
now examined. Sequence x,(n) is zero padded as described in 
Equation 21. A frequency domain series is created using the 
FFT as described by Equation 18. The autocorrelation estimate 
is derived using Equation 19 and is plotted in Figure 13. The 
results of frequency domain based autocorrelation are very 
Similar to the results of time domain based autocorrelation 
using shifted Rayleigh noise. 
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Figure 13. Normalized frequency domain based autocorrelation. 


3. Spectral Interpolation 

With the incorporation of the FFT in modern radar 
receivers, it is possible that the received pulse burst has 
been transformed into frequency data as a byproduct of 
detection. Processing time would be lost in converting this 
data back into a time series to be correlated. Of course, 
there is no guarantee that the dc component has been removed. 
In this section an algorithm is given to perform frequency 
domain based correlation on a Rayleigh distributed noise 
sequence while avoiding the triangular correlation bias 
created by the (possible) dc component. 

A Rayleigh distributed noise series x(n), N points 
long, is transformed into an ordered sequence of N complex 


coefficients using the FFT algorithm. In general, a sequence 
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X(k) of coefficients is indexed from zero to N-1 and is 


described by 


X(k) FFT { x(n) }, (an 


where k = frequency index 


(31) 


The FFT creates the zero indexed term X(0) by summing all the 
elements in x(n). X(0) is real and can be viewed as the dc 
component of x(n). By removing the X(0) coefficient from the 
frequency sequence we have in effect removed the dc bias from 
the time series x(n). At this point, the modified Fourier 
transform X(k) approximates the DFT of a shifted Rayleigh 
series. 

Time domain correlation of an N point series produces 
2N data points. Frequency domain based correlation of an N 
point sequence produces N data points. To force frequency 
domain based correlation to equal time domain based 
correlation (i.e., avoid circular correlation), the number of 
terms in the frequency series must be doubled. The frequency 
series is expanded by storing each pair of complex frequency 
coefficients at a location twice the value of the original 
index. 

As a first approximation, the data points between 
adjacent coefficients are linearly interpolated in the 


expanded array. Figure 14 illustrates the effect of FFT 
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interpolation and the zeroing of the dc terms. Correlation 


results as described by Equation 19 are shown in Figure 15. 
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Figure 14. FFT interpolation, and zeroing of the dc 
coefficient. 


normalized iFFT 2N interpolated autocorrelation 





Figure 15. Normalized frequency autocorrelation through 
interpolation. 
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Clearly, the triangular bias in the correlation output has 
been removed. A small impulse at a time lag of -N is observed, 
which can be disregarded. 

For block processing, typically only correlation 
products within +10% of the zero time lag position are 
considered to be statistically reliable. Therefore, the TDOA 
estimate is formed within this correlation band. The impulse 
at -N falls outside of the TDOA evaluation area and will not 


cause a false detection. 


C. CORRELATION OF PULSES IN ADDITIVE RAYLEIGH NOISE 

Time domain based correlation of a pulse train in additive 
Rayleigh noise should produce the same result as frequency 
domain based interpolation followed by an equivalent 
correlation. A comparison was performed using a 50% duty cycle 
pulse train in additive shifted Rayleigh noise. 

Time correlation was performed on the pulse burst. The dc 
component of the noisy pulse train is removed forming a zero 
mean signal. The time domain based autocorrelation of this 
sequence is shown in Figure 16(a). 

The pulse burst was also transformed into the frequency 
domain using the FFT with the dc component intact. N Fourier 
coefficients are expanded and the X(0), X(+1) coefficients 
zeroed out. The sequence is correlated after interpolation 


using Equation 19 and plotted in Figure 16(b). 
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Figure 16. (a) Normalized time domain based autocorrelation. 
(b) Normalized frequency domain based interpolation and 
autocorrelation. Zero time lag between each channel. 


Clearly, the two correlation algorithms produce similar 
outputs. Both have a dominant peak in their output for TDOA 
estimation. Frequency domain based correlation displays a 
exponential decay of the sidelobes about the main peak. This 
effect is not observed in time autocorrelation, and is caused 
by the first order interpolation used in the frequency domain 
based method. A more robust interpolation algorithm will 
eliminate the exponential decay, yielding the desired linear 
decay. 

Next, time and frequency domain based crosscorrelation are 
performed on two pulse sequences. One pulse sequence lags the 
other by 20 units. Both time and frequency domain based 


crosscorrelation of the two sequences are plotted in Figure 
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17. Clearly, both correlation algorithms produce a peak at 20 
units. The nonlinear decay of the sidelobes surrounding the 


dominant peak is again seen in frequency domain correlation. 
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Figure 17. (a) Normalized time domain based crosscorrelation. 
(b) Normalized frequency domain based interpolation and 
crosscorrelation. 20 unit time lag between each channel. 


In summary, the frequency domain based interpolation and 
crosscorrelation algorithm produces an output comparable to 
time domain based crosscorrelation. It offers a method to 


correlate pulse bursts collected in the frequency domain. 
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V. THRESHOLD DESIGN FOR ENVELOPE CORRELATION 


A. INTRODUCTION 

A threshold at the output of a correlator is used to 
estimate the location of a TDOA peak in a constant false alarm 
rate (CFAR) receiver. The value of the threshold is a function 


of the correlation output mean and variance. 


B. CORRELATION OUTPUT AND CONSTANT VARIANCE 

If the variance of the correlator, when driven by two 
signals at the input, maintains a fixed value about some mean, 
a constant threshold above the mean can be used. This allows 
the design of a constant false alarm rate detector. The 
largest correlation peak that also crosses the threshold is 
defined as the TDOA estimate. The constant threshold design is 


the easiest to implement. 


C. CORRELATION OUTPUT AND TRIANGULARLY SHAPED VARIANCE 

The crosscorrelation of two zero mean noise sequences with 
the same number of data points, produces an output whose 
variance has triangular form. Hence, a constant detection 
threshold is not easily implemented. Instead, a threshold must 
be designed that has a slope that matches the slope of the 


correlation output. The design of a sloping threshold is 
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difficult, particularly if the slope of the correlation peak 
changes with time (i.e., changing noise environment). 

A solution is to mathematically force the triangular 
correlation variance to be flat. This is done by using a 
bowtie correction of the form 1/(N-|0|)?/?. A bowtie correction 
normalizes the correlation output to have a constant (lag 
independent) variance. This design suffers from time changes 
in the correlation output forcing the bowtie correction to be 
recalculated. Note that the amplitude of the TDOA peak will 
itself be dependant on the type of correction to the 


correlation cut. 


D. THRESHOLD BASED ON ZERO LAG 

The zero lag threshold method uses the variance of the 
zero lag correlation output to calculate a constant detection 
threshold. This is indicated in Figure 18. 

The correlation output becomes statistically less reliable 
the further removed the correlation products are from the zero 
time lag position. To use correlation intelligently, the 
correlation output is typically examined +10% of the distance 
from the zero time lag position. Consequently, only 
correlation output +0.1N from zero time lag will be used in 
TDOA estimate simulations. 

The detection method in this thesis uses a constant 


threshold against the correlation output. The fixed threshold 
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Figure 18. Correlation output with zero lag threshold. 


was selected because the correlation output variance is, to a 
first approximation, nominally flat within the £0.1N boundary. 

The detection scenario that drives this thesis defines a 
minimum radar pulse width of one microsecond (usec) and a 
burst length of four millisecond (msec). A 50% duty cycle 
radar pulse is assumed. There are 2000 pulses per burst, which 
when correlated, will generate an output that is 8000pnsec 
(z4000unsec) in length. If ¿10% of the correlation output 
around zero time lag are compared to the threshold value, 


+400U4sec will be tested. For radar pulses traveling at the 
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Speed of light, a maximum distance of +120km is therefore 
examined to form a TDOA estimate. 

Correlation values at lags -0.1N to 0.1N will be compared 
in magnitude to the zero time lag based calculated threshold 


value. 


E. THRESHOLD CALCULATION FOR CFAR 

The detection threshold is calculated using the output 
terms of the crosscorrelation algorithm when noise only is 
injected into the correlator. For convenience, Equation 8 is 


repeated here, 


N-|1|-1 


Lu) c DE. 


I1 "9 


(32) 


x, and y,,, are zero mean Rayleigh (shifted) 
distributed noise sequences 


The correlation function at all lags of interest is thought to 
be Gaussian distributed. According to Equation 32, every 
estimate is the sum of a fairly large number of products 
(i.e., central limit theorem). The average value of the 


correlation output between the +0.1N endpoints is 


*a 


= 1 
t >a + 1 DE PEEL) 


l=-a 


(33) 


where a = number of data points between the 


O and 0.1N lag positions in the 
correlation output 


The variance of r.(1) 1s caleulacce@raen 


42 


2 ah ibe 2 
A y (ero eia Ep) ae (34) 


i=-a 


For Gaussian distributed noise the probability of false alarm 
(Pg) is a function of the detection threshold 


Te 


a a 20% 
Pra 7 ———| e"^* dx 
ü c (35) 





where T - detection threshold 


For a given P,, the threshold can be calculated from 


U= erro O 


(36) 


where erfc = complementary error function 


Equation 36 relates the threshold to the Gaussian distributed 


noise variance. 


F. VERIFICATION OF CFAR THRESHOLD AND ROC CURVE 

Using the above equations, thresholds were calculated for 
the time correlation detector given in Equation 32 for three 
different P,'s. A MATLAB simulation was performed to compare 
the measured false alarm rate of the time correlation detector 
to the theoretical false alarm rate. 

1. P, Simulation 

Two 100 point sequences of zero mean Rayleigh noise 

were injected into the correlation detector. The output of the 
correlation detector between the -0.1N and +0.1N lag points 
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was compared to the calculated threshold to measure the false 
alarm rate. The number of points in the correlation detector 
output that exceeded the threshold for a given P, were 
recorded. The experimental P,, is given as the ratio of 
improper threshold crossings to the total number of monitored 
range cells of the correlator output (see Table 1). For each 
P,, 500 realizations were performed. Both theoretical and 


experimental results are listed. 


TABLE 1 


THEORETICAL PROBABILITY OF FALSE ALARM 


0709/99 OOOO 





The data in Table 1 show that a detection threshold 
can be established for a correlation detector, using the 
assumption of a Gaussian distributed output for P,,’s of 0.01 
and 0.001. For the lower P,, value of 0.0001, the threshold 
calculation fails, but in a positive sense. The measured P,, 
of 0 is below the theoretical P, of 0.0001. Either the 
correlation variance is smaller than Equation 36 predicts, or 


more terms must be summed in the correlation output to better 
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approximate the Gaussian pdf as indicated by the central limit 
theorem. 
2.  TDOA Simulation 

The output of the envelope detector will either be 
Rayleigh or Rician (non-central Rayleigh) distributed. In 
those portions of the data where there is no signal energy, 
the output will be Rayleigh distributed. Where the signal is 
present, the envelope detector output will be Rician 
distributed. 

The purpose of the time domain based correlation 
detector is to measure the TDOA of pulsed sequences embedded 
in Rayleigh noise. Noise in low SNR environments will mask the 
TDOA peak. A MATLAB simulation was performed to examine the 
performance of the correlation detector in a decreasing SNR 
environment. 

Two pulse bursts embedded in additive Rayleigh noise 
were created. The second pulse burst lagged the first by three 
time units. Each pulse burst is a sequence of ten pulses, each 
having a period of ten units. Each pulse in the burst had a 
50% duty cycle. 

The SNR established the amplitude relationship between 
the noise and signal sequences. First, two sequences of zero 
mean unit variance Gaussian distributed noise were created. 
Rayleigh noise was generated from the two Gaussian sequences. 
The Rayleigh noise had a mean of 1.2533 and a variance of 
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0.4292. The theoretical second moment of the Rayleigh noise is 
two and was verified in the MATLAB simulation. The second 
moment is defined as the average power of the Rayleigh noise. 
The amplitude of the nonzero elements in each pulse burst was 
then scaled to meet the desired SNR, using the measured 
average power of the Rayleigh noise. 

A detection threshold was calculated based on shifted 
Rayleigh noise using Equation 36 where 0? was assumed to be 
known (158 Simulation information) In practical 
applications, 0? would be estimated from the data. The SNR of 
the two bursts were varied from OdB to 18dB in the different 
simulations. For each SNR, the two sequences were correlated 
using Equation 32. 

Graphs of the correlation output for a SNR of 1dB are 
given in Appendix C for P,'s of 0.01, 0.001790. 000i 91 
0.00001. The threshold, which is a function of the P,, and 
noise variance, is shown to increase with decreasing Pp. 

The correlation peak at a time lag of three was the 
correct TDOA estimate. If the maximum correlation peak in the 
output was the third lag point, a correct TDOA estimate was 
obtained. For each SNR, 30 TDOA estimates were made. The 
correlation function output 100 data (lag) points per TDOA 
estimate. The TDOA estimate was made on +10 data points on 
each side of the zero time lag position in the correlation 
output. The percent correct TDOA estimate for each SNR was 
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calculated and is presented as a receiver operating 
characteristic (ROC) curve in Figure 19. 

The ROC curve shows that as the SNR increases, the 
performance of the correlation detector increases. The 
performance of the correlation detector in Figure 19 can be 
improved if multiple sequential looks at the same SNR were 


allowed. 
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Percentage of estimated TDOA correct versus 


Figure 19. 


domain based correlation detection. 


SNR for time 
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VI. FFT ERROR ANALYSIS 


A. INTRODUCTION 

This chapter discusses several sources of processing 
errors occurring within the FFT, and the effect of these 
errors on the FFT output signal-to-noise ratio. 

The computation of a FFT using fixed point arithmetic, 
generates internal processing errors. These errors are due to 
the fixed number of bits in the trigonometric look up table, 
the truncation of arithmetic results during addition and the 
scaling of results during the individual butterfly 
calculations. 

The error sources can be modeled as independent additive 
noise components at the output of the FFT. Their net effect is 
to reduce the ideal SNR at the output of the FFT. 

Each of these processing error sources will be described 


below. The material in this chapter borrows from [Ref. 7]. 


B. COSINE/SINE TABLE NOISE 
The computation of a butterfly algorithm in the FFT 
involves multiplication by the complex coefficient (twiddle 


factor) 
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cos (^m) - j sin( m) (37) 
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A FFT butterfly is shown in Figure 20. 





Figure 20. FFT butterfly calculation. 


These trigonometric coefficients are precalculated, quantized 
to (B) bits (where B is the size of the computer word 
including sign bit), and stored in memory for future table 
lookup during the FFT processing. The coefficients W° = 1 and 
Ww‘ = -j can be obtained in an exact form and therefore do not 
contribute to the quantization noise. Quantization of the 


other (twiddle) coefficients stored in the sine and cosine 
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table introduce noise into the output of the FFT. The noise 
power created by truncation of the sine and  cosine 


coefficients is defined by [Ref. 7] to be 


2 le 
o LE x 
(38) 
where 0%, = trigonometric noise power 


ty 
i 


number bits in computer word 


E, TRUNCATION NOISE 
The DFT of a finite duration sequence í(x(n)) is defined in 


Equation 18 and is repeated here for convenience 


N-1 2R 


X(k) = y Me an O <s kK<N-1 . (39) 


n=0 


The product formed requires four real multiplications because 
both the exponential and the input data are complex numbers, 
B bits in length. Each multiplication can produce a result of 
2B bits which must be truncated from 2B to B bits, hence there 
are four quantization errors for each complex valued 
multiplication. The four quantization errors can be described 
as four independent noise sources. The truncation noise power 


is defined by [Ref. 7] as 
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o? - [osos on 
truncation noise power 
number of bits in product after truncation 


where os 
B 


D. SCALING NOISE 

The FFT as an algorithm processes a data vector of length 
N by successive passes at the vector. During each pass, the 
algorithm performs N/2 butterfly operations. Each butterfly 
retrieves two complex numbers from memory, performs the 
butterfly computation and returns two complex numbers to the 
same memory address. The complex numbers returned to memory by 
the butterfly can have a greater magnitude than the two 
complex numbers initially fetched from memory by the 
butterfly. In order for the results of the butterfly operation 
to fit in the fixed word length of the memory of the computer, 
the results must be truncated (scaled). Scaling is performed 
by dividing by two the entire data array, and is implemented 
in software by a right shift and an increment of the arrays 
exponent register. 

During the (m+1) pass of the data, the butterfly algorithm 
selects two data points A(m) and B(m) and returns to memory 
A(m+1) and B(m+1). The truncation error results from the 
addition of two numbers of like sign in the upper part of the 
butterfly algorithm or the subtraction of two numbers of 
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opposite sign in the lower part of the butterfly algorithm. 
The input to the butterfly is fixed to have maximum real and 
imaginary values of +1. The output of the butterfly is fixed 
to have a maximum value of +2, double the input value. To 
prevent overflow of the butterfly output, either a prescale by 
1/2 is required prior to entering the butterfly or an extra 
bit must be available in memory to accommodate the possible 
gain of two. 

One method, automatic prescaling, truncates the least 
significant bit of the butterfly output and represents 
processing noise. In the case for which no scaling is 
necessary for the pass, this truncation represents significant 
processing error (noise). 

A second method, data dependent scaling, is performed only 
if a butterfly in the data pass overflows without scaling. If 
any word in memory has an overflow bit set, all words are 
shifted right as they are fetched from memory for the next 
pass of the FFT. The total number of right shifts, p, executed 
during the FFT process is used at the output of the FFT as a 
Scale factor of 2P. This scale factor is used at the output of 


the FFT so that the total processing gain is maintained. 
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The noise power created by the scaling of the data is 
defined by (Ref. 7) 
2. 1, 1, 2e l >-2(B-1) 
E + =2 
LE 4 


(41) 
where 0% = scaling noise power 


E. OUTPUT SIGNAL TO NOISE RATIO 


Figure 21 shows an error model of a butterfly in an FFT. 





Figure 21. FFT error model. 
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The input data terms to a butterfly can be described by 


A(m) = A(m) + N,(m) 
B(m) = B(m) + N,(m) 
ME : (42) 
wnere N,, N, = addi€ive nolse terms 
A(m), B(m) = input data terms 
A(m), B(m) = noise corrupted data 


Without the effects of noise the butterfly algorithm produces 


two outputs described by 


A(m+1) = A(m) + W(m) B(m) 
B(m+1) = A(m) - W(m) B(m) (43) 
where W(m) = twiddle factor 


The error due to the A component will be calculated. The error 
due to the B component is statistically equivalent to the 
error in the A component. Including the effects of noise, the 


output of the butterfly can be described by 


A(m+1) = [A(m) + N,(m) + q,(m)] + [B(m) + N,(m) +q,(m)]' 
[W(m) + q,(m)] + g,(m) 
where q,(m) = scale noise 
q,(m) = sine/cosine table noise 
q,(m) = truncation noise (in multiply) 


(44) 


The total power at the output of the butterfly is then 
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E{A(m+1) Ä’(m+1)} = E{A(m)A’(m)} + E{B(m) B*(m)} + 
E{N,(m)Na(m)} + E{N,(m) Ng (m) } 


: A 45 
Etqs(m qi(m) + Etgstmajimy + 9) 
E(B(m) B' (m) d] Eig, (na. (me 
E{g-(m) gr (m) } 
Redefining the terms of the total power leads to 
E{A(m+1)A*(m+1)} = S2(m) + S?(m) + Ng(m) + Ng(m) + 
202(m) * S?(m o,(m) * o£(m) 
where S?(m) = input signal power (46) 
o2(m) = scaling noise power 
o%(m) = truncation noise power 
o%(m) = sine/cosine table noise power 
Combining like terms results in 
E[A(m+1) A" (m+1)) = 25?*(m) + 2N¿(m) + 205(m) + (47) 


s?(m)oy(m) + o7(m). 


This equation describes the total power as a function of the 
signal power and the individual independent uncorrelated noise 
terms. The total power output from the butterfly can also be 
described as a function of its input signal power and the 


total additive noise power generated within the butterfly as 


E{A(m+1)A*(m+1)} = E{A(m+1)A*(m+1)} + (48) 
E(N¿(m+1) Na (m+1)) 


Redefining these terms 
E{A(m+1)A*(m+1)} = S?(m+1) + Ni(m+1). (49) 
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Equation 47 and 49 are equivalent. The signal and noise power 
for the (mtl) pass can be formed from those terms in the 


proceeding pass 


S*(m+1) + Ni(m+1) = 25?(m) + 2N)(m) + 202(m) + (50) 
S?(m)o;(m) - o£(m). 
Equating the signal terms in Equation 50 yields 
Sm 2S A (51) 
Equating the noise terms in Equation 50 yields 
Ni(m+1) = 2Ni(m) + 202(m) + S?(mo%(m) * o?(m). (52) 


Equation 51 can be rewritten to show how the input signal 
power increases as a function of m as it passes through the m'^ 


Stage of a butterfly of an FFT 


See aus 


(53) 
here we modeled S?(0) = S?. 
The signal power input to the butterfly is defined as 
Sg? = 15-2 -b . (54) 
9 


Rewriting Equation 52 to include this exponential increase of 


S? (m) gives 
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Ni(m*«1) -2Nz(m) * 20$(m) * 27S?gí(m) * 0%(m). (55) 
Equation 55 shows the relationship between the noise power of 
a butterfly stage as a function of the noise power from the 
proceeding butterfly stage. The noise and signal power can now 
be calculated using Equations 53 and 55 for m passes through 
a butterfly. Once calculated, the signal to noise ratio can be 


formed. 


F. RIGHT JUSTIFIED DATA 

The length of a word in the computer is defined to be 16 
bits (B) and accepts a 12 bit (b) data word. Right justified 
data places the least significant bit (LSB) of the b bit data 
word in the LSB of the B bit processor word. The output noise 
to signal ratio is now calculated for ten passes through a 
butterfly using right justified data. The following 
assumptions are used in the calculation; 

1. The input signal is incoherent. 

2. Scaling is performed every other pass after the most 
Significant bit (MSB) of data reaches the next to MSB 
of the processor. 

3. The first two passes do not use the trigonometric 
table. 

After the tenth pass (i.e., FFT size of 1024), it is shown 


(Appendix D) that the noise to signal power is 
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N? N: On c? 
Eu MM LM EM ee (56) 
3 4 


B and b are previously assumed to be 16 and 12 bits in length 


respectively. Substituting these values into Equation 56 





Y; s: ae Ead 
SONNO ) 2 ? £4 55022 EE np . (57) 
S? 13-8 5 4 13-8 a 13-8 
9 9 9 
The noise to signal ratio reduces to 
N: 
<a (10) = -67.5 db - 95.1 dB -71.1 dB -79.1 dB (58) 
SPUC crag truncation scaling 
noise noise noise noise 


Clearly, the dominant degradation noise is the truncation 
table noise (-71.1dB), followed by the scaling noise to signal 
ratio (-79.1dB). For right justified data, the signal to noise 


ratio is calculated to be 65.71dB. 


G. LEFT JUSTIFIED DATA 

Left justified data places the most significant bit (MSB) 
of the b bit data word in the MSB of the B bit FFT processor 
word. The output noise to signal ratio is now calculated for 
10 passes through a butterfly for left justified data. The 


following assumptions are made 
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1. The signal is incoherent. 
2. Data is left justified on input and scaling is 


performed every other pass starting with the first 


pass. 
3. The first two data passes do not use the 
trigonometric table. 


It can be shown that after the tenth pass through the 


butterfly that the noise to signal ratio is 


2 2 2 
E NOE = E + 40% + 1 9T + 7 (59) 
S? Sg? v 4 Sg? 


B and b were previously assumed to be 16 and 12 bits 


respectively. Substituting these values into Equation 59 we 


obtain 





2 i278 1l 7 5-32 2 [27 
Aa. - n EM > - (60) 
2 T Tum EE 
=? — =2 
9 9 9 
The noise to signal ratio reduces to 
Na 
cL = -67.5 dB -95.1 dB -89 MDE UN (61) 
- input trig truncation scaling 
noise noise noise noise . 


Clearly, the worst noise contribution to the SNR is the 
scaling noise (-77.1dB), followed by the truncation noise (- 
89.1dB). For left justified data, the signal to noise ratio is 


calculated to be 67.0dB. 
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Further analysis shows that for large FFT's (i.e., in the 
range of 2 and above), the trigonometric table noise is the 
dominant noise source. A simulation was performed using an 
AMDAHL machine to examine the error generated in performing a 


FFT and an iFFT. Figure 22 shows the experimental setup. 


AMDAHL Machine 


Y 
save error in table 





Figure 22. FFT processing error flowchart [Ref. 8]. 


The error that was generated was tabulated. Figure 23 


illustrates the error as a function of FFT length. 


H. SUMMARY 

For either left or right justified data, the dominant 
noise at the output of the 1024 point FFT is the scaling or 
truncation noise. The SNR at the output of the FFT, for left 
justified data is 67.0dB and exceeds the SNR for right 
justified data which is 65.71dB. Clearly, left justified data 
maximizes the computational SNR at the output of the FFT for 
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MAXIMUM ERROR IN SINGLE 
PRECISION FFT 


A 
10C,000—— 
di / ` 
| / 
CVV a 
C Ne 
(C "un 
SO 2000088 — 
cc | 
tr 33.000 s 
W c 00 





5 6 7 8 9 10 1 12 13° 14 3516 LZ 


FFT LENGTH 


Figure 23. FFT error as a function of transform length 
(Ref. 8]. 
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the values of b, B and N used in the above examples. 
For a right justified data set transform size of N-2? 
(8192 points), the noise to signal ratio after the 13'^ pass 


can be written 


2 


Ma) - 2 nó + (225 +20 + 212) 5?%07 + 20470%+ 17004 (62) 
This equation reduces to 
2 2 2 2 
Na (43) ME 00.02, al NEN US (63) 
S? g? 4s? 48.2 g? 


Comparing Equation 56 to Equation 63, we see that by 
increasing the transform size from 2! to 2? increases the 
noise to signal ratio. The trigonometric table noise is 
increasing faster than either the scaling or truncation noise. 
The signal to noise ratio is reduced, from 65.71dB to 65.68dB. 
As the transform size is increased to accommodate larger data 
sets, the trigonometric table noise becomes the dominant error 


term. 
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VII. SUMMARY 


A. CONCLUSIONS 

A frequency domain based algorithm was developed and 
tested to estimate the differential arrival time of a pulsed 
radar signal collected by two passive sensors. For 
convenience, the actual implementation was performed through 
time domain processing even though frequency domain processing 
is advocated. The performance of the algorithm is 
characterized by a ROC curve as a function of SNR. For 3dB SNR 
and a P, of 0.01, the probability of making a correct TDOA 
estimate exceeds that of an incorrect estimate for a single 
look. At 18dB SNR, the probability of making a correct TDOA 
estimate is 100 percent for all P,,’s. If multiple looks (i.e., 
multiple bursts) are allowed, the probability of making a 
correct TDOA decision at each SNR will increase. 

An I/Q demodulator is assumed in the radar receiver. The 
pdf of the signal driving the correlation detector is 
determined from where it is obtained in the I/Q receiver. The 
pdf, depending on selection, will be zero mean Gaussian versus 
non-zero mean Gaussian, or chi-squared versus non-central chi- 
squared, or Rayleigh versus Rician (non-central Rayleigh). 
This thesis assumes an envelope detector at the output of the 
I/Q demodulator. The envelope detector output has a Rayleigh 
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or Rician pdf and drives the correlation detector. The 
calculation of a CFAR threshold assumes a Gaussian pdf at the 
output of the correlation detector. For small time lags, 
summation of the output of the correlator may produce a 
probability distribution that deviates from Gaussian. This 
deviation would produce a biased threshold calculation. To 
minimize the bias, a sufficient number of terms must be summed 
at the output of the correlator to allow the Gaussian 
approximation. 

If the received pulses are collected in the frequency 
domain, a spectral domain based correlation algorithm can be 
implemented. An algorithm is given in this thesis. 

For any digital signal processing algorithm that uses the 
FFT, processing errors must be considered. For large transform 
sizes the trigonometric noise power dominates. The length of 
the trigonometric coefficient word affects the degradation of 
the SNR at the output of the FFT. The larger this word the 
smaller the noise power. For a given transform size, left 
justified data will have a higher output SNR than right 
justified data. 

For a correct TDOA estimate, the correlation algorithm 
requires the reception of the leading pulses in the radar 
pulse burst. If the pulses are received in an adequate SNR, 
the correlation algorithm will produce a well defined peak 
allowing the TDOA estimate. Should the environment degrade the 
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received signal (i.e., destructive interference due to 
multipathing, weather or terrain), the correlation algorithm 
should be discarded in favor of the traditional angle of 
arrival (AOA) algorithm (1.e., maximize received energy). 

If the radar uses a staggered PRF, the reception of the 
first pulse or several pulses has little impact on the 
position of the correlation peak. Under this condition, the 


TDOA algorithm is superior to the AOA algorithm. 


B. RECOMMENDATIONS 

A ROC curve should be constructed for the correlation 
algorithm used in this thesis for multiple TDOA looks. An 
intensity display could be designed that would plot as a 
function of time (i.e., snapshots), those lag points chosen as 
TDOA estimates. Patterns displayed on the intensity plots 
would allow the user to visually determine the correct TDOA 
(i.e., incoherent averaging). 

The performance of the spectral domain based correlation 
algorithm developed in this thesis must be further quantified. 
Sets of ROC curves should be obtained for different SNR’s 
(1.e., SNR variations between channels). 

The simulated FFT error graph should be validated by 
measuring the error performance of a dedicated properly 


dimensioned FFT. 
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APPENDIX A. CORRELATION MEAN AND VARIANCE 


Two zero mean, independent time series are the inputs to 
a correlator. This appendix will derive the expected value and 
variance of the crosscorrelation function output. In this 
thesis, the noise is zero mean (i.e., shifted) Rayleigh noise. 
The crosscorrelation function r,,(1) is defined in Equation 


8 and is repeated here for convenience 


N-1-|1| 
ZEN rev © (64) 


1=0 


A. EXPECTED VALUE 
x, and y, are assumed to be independent, zero mean 


sequences 
E{x,} = Ely) = 0, (65) 
The expected value of r,,(1) is 


N-1-|1| N-1-|1| 
E(r,()) 9 ECUCO, xyuüib- Y Elixy.a 
1=0 1=0 


N-1-|1| 66 
Y EG El.) E 
1=0 
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B.  VARIANCE 


The random variable z has a variance defined to be 


2 


o; = El(z Eis 
because E{z} = 0 for a zero mean sequence, (67) 
t = ez 


The variance of the crosscorrelation function r,,(1) is then 


3 N-1-{1| N-1-|1| N-1-11| 
o CI ELE xy) = El py "Ya UNE > X1Y1«1 = Xi: Via 
i i p (68) 


N 


N-1-|1| N-1-11| 


Eid AA A 
izo 3% 


xx, and y, Yyı are two groups that are independent of each 
other. The expectation operation can be applied to each group 


individually. 


N-1-|1| N-1-|1| 


Ort) = 2 2 Ex x, Evo M (69) 
=0 =0 


The two indices i and j define a matrix of values. Two 
contributions have to be considered. 

Contribution 1. When i-j the summation occurs along the 
diagonal of the matrix. Only one summation is used and 
essentially the terms are being squared. These squared terms 
are not independent. To simplify the computation define m 


equal to i and j 
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N-1-111 
NS E ME [Y (70) 


m=0 
Because x, and y, are zero mean sequences, 
El) O SU > 
and 0% = 0, = 0? 
Substituting 0° into Equation 70 


N-1-|2| 


N-1-|1| 
yoo; = of Yt 
m=0 m=0 


Cy 


(71) 


o*{ (N-1-|1]) +1] 


D 


Contribution 2. When izj all the terms except those on the 


diagonal are summed. These terms are independent. For izj 


N-1-|1| N-1-12|] 


C, = EY E.) 
CM A ge (72) 


0, using Equation 65 


Therefore, Equation 69 becomes O0; = C, + C, (73) 


of(N - | 1|) 
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C. CROSSCORRELATION STATISTICS 

Each data point output by the correlation functo 
represents a summation of terms. Using the central limit 
theorem, the output of the correlation function is assumed to 


have a Gaussian distribution. 


In conclusion (OD NO OA ra 


where N = Normal (Gaussian) distribution 
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APPENDIX B. MATLAB CODE 


A. USERS GUIDE 

The frequency domain based correlation software given in 
this appendix is written in Pro-Matlab (version 3.5h). All 
Simulations were run on a Naval Postgraduate School Sun work 


Station using a UNIX operating system. 


B. CORRELATION MATLAB CODE 


EEEEEEEESES EES EERE EE EEE EE EEE EE ES EEE EEE EEE ES ES EE EEE ESE ESES 
$ Program Name : FreqCorr7.m 12 May 91 
$ 


ES EEEEEEEEEESERE EEE EE EEE EE SES EE EEL EE EEE EEL ER EEE EEEEE EEE ESEEES 
clear 


clg 

CEC 
EEEEELEFESEEEEEEE EEE EEE EEE EEE EEE EEE EEE ES EEE EEE EEE EEE EEE EEEES 
$ 

% PULSE CORRELATION 

$ Version 7.0 

$ 

% 1. 1 Channel vs Reference Channel 

$ or 

$ 1 Channel vs Second Channel. 

$ 

$ 2. Rician Signal and Rayleigh Noise pdf's 

$ (Low SNR signals have power approximately 2.) 

$ 

$ 3. 0 dc Correlation. Subtract dc from both signals. 

$ 

$ 4. Normalize correlation output in Fourier domain. 
EEEEEEEESESEEEEEEEESEEEEEF EES ESEFESE EEE EEE EES EL ES EL EEEEEES 
% This MatLab program will either 

% 1. Correlate a received pulse sequence against a reference 
$ pulse sequence. The reference sequence parameters are 

$ specified by the user and are used to construct the 

$ sequence. 

$ or 

$ 
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% 2. Correlate two received pulse sequences against each 

$ other. 

$ 

$ 

$ Correlation of the pulse sequences is accomplished by 

$  multipling in the frequency domain one spectrum against 
% the conjugate of the other spectrum. An inverse FFT is 

% performed on the result to yield the time difference 

$  (TDOA) between the two sequences. 
EEEEEEEEEEEEEEEEEEE EEL EEE E EEE EEE EEE EEE E EEE EEE EEE EEE EEE EEEES 
CIC 

echo on 

EEEEEEEEEEEE EEE EEE EEE EE ESE EEE EEE EEE EERE EEE EEE EEE EEE EEE EEEES 
% SELECT CORRELATION OPTION 

% 

% Select either 1 or 2 

% 

$ 1. Time delay in one receiver channel measured against a 
% reference signal. 

% 

% 2. Differential time delay between two channels. 

% 


EEEEEEELEE EEL EE EEE EEE EEE EEE EEEEELEEEEEELEEEEEEEEEEEEEEEEEEES 
echo off 


option = input(’Select correlation option 1 or 2 ; 2 
elle 

echo on 
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 
% Option 1 : Time Delay In One Channel 

$ 

$ Reference Pulse Train Parameter Selection. 


ESEEEEEEEEEEEEEEEEE EERE EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 
echo off 
if option== 


wi-input( 'Pulse Width (seconds) : I os 
pr-input( 'Pulse Period (T seconds) : o 
nuzinput( 'Number of pulses in pulse train : $ 
NumberPointsRef=pr*nu; 

end 


eZ 


ole 

echo on 

EFESEEEEEEEEE EE EEE EE EEE EE EEE EEE EEE EEE EEE EEEEEEEEEEEEEEEEEEES 
% 

$ 

$ Delayed Pulse Train Parameter Selection 

% 

% 
EEEEEEESEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 
echo off 


dly=input (’Channel 1 delay tau (sec.) of pulse train : ny? 
wi2-input('Pulse Width (seconds) : adc: 

pr2-input('Pulse Period (T seconds) : zy 

nu2-input (/Number of pulses in Delayed Pulse Train : E 


NumberPointsDel=(pr2*nu2)+dly;  $calc. nmbr pts sig + delay 


ee 

if option== 

dlyCh2-2input (/Channel 2 delay tau (sec.) pulse train : 25 
wiCh2=input (“Pulse Width (seconds) : os 
prCh2-input('Pulse Period (T seconds) : 2 

nuCh2=input (‘Number of pulses in Delayed Pulse Train : B 


NumberPointsDel2= (prCh2*nuCh2) +dlyCh2; $nmbr pts sig + delay 
end 


$1533151515359595593959595939595955395959595959595959959559955939959959939959595939595939959595593** 
Zero pad both sequences to N3=N2+N1-1. Must then make N3 
meet the simple formula N3 = 2%m to allow speedy 
computation of the FFT. 
WARNING : IBM 80286 MATLAB WILL NOT ALLOW VECTORS GREATER 
THAN 4000. So an input pulse train of 128 will exceed the 
System. Solution is to use the SUN workstation or 386/486. 
$5$5525255159554355439599539595395954595995959595959595959595959595959595954595959595159595959595*559593595315 
if option==1 
ZeroPadPoints=NumberPointsRef+NumberPointsDel-1;%N3=N2+Nl-1 
for m=2:1:19; $2” (19) = 524288 point FFT max 
if 2^m >= ZeroPadPoints, 
ZeroPadPoints-2^m; 
break; $when find a 2^m power, exit loop 


oe ge ge oe ge age 


end 
end 
end 
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if option== 
ZeroPadPoints=NumberPointsDel2+NumberPointsDel-1;%N3=N2+N1-1 
for mz2z:1:0195 $2^(19) = 524288 point FFT max 
if 2”m >= ZeroPadPoints, 
ZeroPadPoints-2^m; 
break; $when find a 2^m power, exit loop 
end 
end 
end 


ESSEEESEEEESEESEESESE Create Reference Pulse Train  $$$$$55555454594*5* 
if option== 
ReferencePulse=zeros (1:ZeroPadPoints); $zero pad past pr*nu 
for PulseCounter-1:pr:NumberPointsRef  $steps of Ref. period 
for CutUpPulse=0:pr; 
if CutUpPulse<=wi 
ReferencePulse (PulseCounter+CutUpPulse)=1.0; 
end 
end 
end 
Refdc=mean (ReferencePulse (1:NumberPointsRef)); 
ReferencePulse (1:NumberPointsRef) =ReferencePulse (1:NumberPoi 
ntsRef)-Refdc; 
end 


%%%3%3%% Create Channel 1 Delayed Pulse Train %%%3333%33%3%% 
DelayedPulseTrain=zeros (1:ZeroPadPoints) ;%0 pad past pr2*nu2 
for PulseCounter=(dly+1) :pr2: (pr2*nu2+(dly+t1)) ;$pr2*nu2=Nmbr 
for CutUpPulse=0:pr2; %...pts in Del. Pulse Train 
if CutUpPulse<=wi2 
DelayedPulseTrain (PulseCounter+CutUpPulse) =1.0; 
end 
end 
end 


if option== 
%$%%%%%E% Create Channel 2 Delayed Pulse Train 33%33%3%3%3%3%%% 
DelayedPulseTrain2-zeros(1:ZeroPadPoints); $0 padpast pr2*nu2 
for PulseCounter=(dlyCh2+1) :prCh2: (prCh2*nuCh2+ (dlyCh2+1) ); 
for CutUpPulse=0:prch2; 
if CutUpPulse<=wiCh2 
DelayedPulseTrain2 (PulseCounter+CutUpPulse)=1.0; 


Do 


end 


end 
end 
end 
Ele 
echo on 
EEEEEEEEEEEEEEFEEEEE EE EE EEE EEE EE EE EE EEE EE EE EEE EEE EEEEEEEEESES 
% Signal to Noise Ratio Calculation 
% Define Rayleigh noise power to equal 2. 
% Rician pdf created for sinusoid + Gaussian noise. 


EEESEEEEEEE EEE EES EEE EE EEE E EES EEE EEE EE EEE EEE EEE EEE EEEEEEEEEES 
echo off 

SNR=input (’ Desired SNR :  '); 

NoisePwr=2; $rayleigh noise pwr =2 


EEESEEEEEESEEEEESEEEEEE EEE EEE EE EEEEEEESEEE EEE EEE EEE EEE EEE EES 

% Sum Channel 1 & 2 signals plus their initial delay. 

EEEFESEEEEEEESEEEESESEEE EE EEE EE EEE EEE EE EEE EEE ES EEE EEEEEEEEES 

SigEnergy=0; initialize Channel 1 energy to 0 

for index=1:NumberPointsDel;%sum over Ch. 1 delayed signal 
SigEnergy=SigEnergy+ ( (DelayedPulseTrain (index) ) .*%2); 

end 


if option== 

SigEnergy2=0; initialize Channel 2 energy to 0 

for index=1:NumberPointsDel2;%sum over Ch. 2 delayed signal 
SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2 (index) ) .%2); 

end 

end 


EEESESESEEEEE EEE EEE ES EEE ES ES ES EE EEE EEEEEEEEEEEE ES EEE EEEESEES 
$ Calc. Ch. 1 signal amplitude to meet required input SNR. 
EEEEEELELEEEEEEEEEEEESEEEE EE EE EE EE EE EEEE EE EE EE EE EE EEEEEEEES 
SNR10=SNR/10; 

SNRLinear-z10^SNR10; $SNR in linear units 
SigAmplitude=sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 
gy); 


if option== 
ESEEEEEEEEEEEEEEEEEEEEEEEEEEEEE EEE EEE EEE EEE EEE EEEEEEEEEEEES 
$ Calc. Ch.2 signal amplitude to meet required input SNR. 
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$9$**5**55959*559555599599959*9*959*959595939955959599*9959595**9*29*9*9999*9*9**9«9€* 
SNR10=SNR/10; 


SNRLinear-10^SNR10; $SNR in linear units 
SigAmplitude2=sqrt (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 
ergy2); 

end 


EFESEFSEREEE EEE EEE EES EEE ESE EL EEE EEE EEE EE EEE EEE EEE EE EEE EEEES 
% Create Rayleigh Noise matrix for Channel 1. 
ESEFEEEEEEELESEEEEEEEE EE EEE EF EEE EEE EEE EE EEEEEEEEEEE EEE EE ES 


rand (' normal’); *Gaussian mean=0, variance=1 
for index-1:NumberPointsDel, 

GaussNoisel (index) =rand; $N (0,1) 

GaussNoise2 (index) =rand; $N(0,1) 
end 


RayleighNoise=sqrt ( (GaussNoisel) .*2+(GaussNoise2) .*2); 


SEEEEEEEEEEEE EE EEE EEE EEE EE EE EEE EEE EE EEE EEE EE EEE EEE EEEEEEEES 
% Create Rayleigh Noise matrix for Channel 2. 

EEEEEEEEEE EEE EEE EEE EEE EE EEE EE EEE EE EEE EEE EE EEE EE EEE EEEEEES ES 
if option== 


rand(’normal’); %*Gaussian mean=0, variance=1 
for index=1:NumberPointsDel2, 

GaussNoise3 (index) =rand; $N(0,1) 

GaussNoise4 (index) =rand; %N (0,1) 
end 


RayleighNoise2=sqrt ( (GaussNoise3) .*2+(GaussNoise4) .%2) ; 
end 


EEEEELEEELESESEE EEE EE EEE EE EE EEE EEEEEEEEEEEEEEEEEEESEEEEEES 

% Make Rician amplitude of Channel 1 delayed pulse train 

% and remove the dc component of entire signal. 

$52522$$59$2551595252552529525515525255315523155525543*52545255155515255155451554355*35593531 

SE=0; 

for index=1:NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain (index) == $if one, make Rician 


DelayedPulseTrain(index)-DelayedPulseTrain(index).*rician(Si 
gAmplitude); 
SE- (DelayedPulseTrain(index)).^2-tSE; 
end $end Rician modification loop 
end 
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SP=SE/NumberPointsDel; 


$%%%%%% Add noise to Channel 1 Delayed Pulse Train %%%%%3%% 
for index=1:NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain(index)==0 %no Rayleigh noise to signal 


DelayedPulseTrain (index) =DelayedPulseTrain (index) tRayleighNo 
ise (index); 

end 

end 


$%%%% Remove dc component of Channel 1 Pulse Train %%%% 
Chldc=mean (DelayedPulseTrain (1:NumberPointsDel)); 
DelayedPulseTrain(1:NumberPointsDel)-DelayedPulseTrain(1:Num 
berPointsDel) -Chldc; 


if option== 

EEEEEEEEEEEEESEEEEEEEEEEEEEESEEEEEEEEEEEEEEEEEEEEEEEEEEEESS 

$ Make Rician amplitude of Channel 2 delayed pulse train 

% and remove dc component of entire signal. 

EEESEEEEEEEEEEEESESEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEES 

SE=0; 

for index-1:NumberPointsDel2, %1 --> delay + delayed signal 
if DelayedPulseTrain2 (index) == $if one, make Rician 


DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) .*rician( 
SigAmplitude2); 
SE= (DelayedPulseTrain2 (index) ) .*2+SE; 
end tend Rician modification loop 
end 
SP=SE/NumberPointsDel2; 


%$%$%%% Add noise to Channel 2 Delayed Pulse Train %%%%%%%% 
for index=1:NumberPointsDel2, %1 --> delay + delayed signal 
if DelayedPulseTrain2 (index) ==0 %no Rayleigh noise to sig 


DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) tRayleigh 
Noise2 (index) ; 

end 

end 


%$%%%% Remove dc component of Channel 2 Pulse Train %%%% 
Ch2dc=mean (DelayedPulseTrain2 (1:NumberPointsDel2)); 
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DelayedPulseTrain2 (1:NumberPointsDel2) =DelayedPulseTrain2 (1: 
NumberPointsDel2) -Ch2dc; 


end 
subplot (221) %two rows, two columns 
time=(0:1:ZeroPadPoints-1); $Start time axis at 0 


1f option==1 

%EEEES$%ESEEEEE Plot Reference Pulse Train  BE¥EEEEEEEEEEEEEES 
plot (time, ReferencePulse); 

title(’Reference Pulse Sequence’); 

xlabel(’time’); 

gua 

end 


%%%% Plot Channel 1 pulse train + Rayleigh/Rician noise%%%% 
topDel=1.1*max (DelayedPulseTrain) ; 

bottomDel=1.1*min (DelayedPulseTrain) ; 

axis([0 NumberPointsDel bottomDel topDel]) 

plot (time, DelayedPulseTrain) ; 

title(’Channel 1 Pulse Sequence + Rician/Rayleigh Noise’); 
xlabel('time'); 

grid 


if option== 

%$%%% Plot Channel 2 pulse train + Rayleigh/Rician noise %%% 
topDel2=1.1*max (DelayedPulseTrain2); %10 percent headroom 
bottomDel2=1.1*min (DelayedPulseTrain2); 

axis([0 NumberPointsDel2 bottomDel2 topDel2]) 

plot (time,DelayedPulseTrain2); 

title('Channel 2 Pulse Sequence + Rician/Rayleigh Noise’); 
xlabel('time'); 

grad 

end 


EEEEEEEEEEEEEEEEEEREEEEEEEEEEEEEEE EEE EEE EEEEEEEEEEEEEEEEEEES 
% Correlation 
ESESEEEEEEEEEEEEEEEEEEEEEEEEEE EEE EEE EEEEEEEEEEEEEEEEEEEEEEES 
if option== 

fR=fft (ReferencePulse) ; 

fD-fft(DelayedPulseTrain); 

end 
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if option>= 

fR=fft (DelayedPulseTrain); 
fD=fft (DelayedPulseTrain2); 
end 


Z=fD. “conj (fR); 


$5$525252555959595995959595959995959595959595959595555555555555*5 
$ Calculate normalizing coefficient 
EEEEEEEEEEEEEEEEEEEE EEE EEE EEE EEEE EEE EEE EEES 
Norml=0; 
Norm2=0; 
for k=1:1:length(fD) 

Norml= (abs (£fD(k)) .”2)+Norml; 

Norm2= (abs (fR(k)).^2) +Norm2; 
end 


Norml=sqrt (Norml); 
Norm2-sqrt (Norm2); 


NormCoeff=Norml*Norm2; 
z-z./NormCoeff; 


z-z*length(fD); 


ifftz=ifft (Zz); 
magz=abs (ifftz); 


%%%%% Flip vector for appearance 3%%%3333%3%%%3%%3% 
Topmagz= (length (magz) /2+1); 
Bottommagz=length (magz) /2; 

magzLen-length(magz);  $measure length of magz 
Transmagz-[magz (Topmagz:magzLen),magz (1:Bottommagz)]; 


%$%%E%%% Make -time to +time axis %%%%%3%%%9%%%%%%% 
minustime=- (length (Transmagz) /2) ; 
posittime= (length (Transmagz) /2)-1; 
timel=(minustime:1:posittime) ; 


topmagz=1.1*max(magz); tadd 10 percent head room 
axis([minustime posittime 0 topmagz]) %scale axis 


subplot (212) 
plot (timel, TransmagZ) ; 
title(’FFT Correlation : Channel 1 vs. Channel 2’); 


T9 


xlabel(’time lag’); 
grid 
gtextí('t7dB SNR) 


lj4print %Spanagel Room 427 Laserjet 


axis([1 2 3 4]); %reset axis scaling 
axis; 


EFEEEEEEEEEEEEE EEE EE EEE EE EEE EEE EE EEE EEE EEE EEEEEEEEEEEEEEEEES 
$ Program Name : CorrPfaLoop.m 19 July 91 
*1$3*551*11*1í$1í$15**5454*93954595959599595959951595955959595959595953595955959595959559595595959595*9 
$ Purpose : To compute the Pfa for Signal / Noise only for 

$ then correlation function. This program creates zero mean 
% RAYLEIGH noise. 
€$5$*5151333559554595954535959595959559595959559595955954954555954555959595959599559595959595*5** 


Clear 

CEG 

subplot (221) 
rand(’normal’) 


PFA-1.0; twill be divided by ten to start 
EE EEEEEEEEEEEEEEEEEEE EEE EEE EEEES 

for PFACOUNTER=1:1:4 

KEE EEEEEEEEEEEEEEEEEEEEEEEEEEEEES 

PFA=PFA/10; &Pfa = 0.1 0.01, 0.001, 09017 


EEEEEGEEEEEEEE EEE EEE EEESEESS 
for SNR=0:1:18 Outer loops through all SNR's. 
$1535551531$1$1515495959595959595959595959559** 


Iterations=30; $30 data points per SNR point 
ESESEEEEEEEESEEEEEEEEEEE EES 

for loopzl:l:Iterations $inner loop, creates groups of +- a 
EEEEEEE EEE EEE EEE EE EEE EEEEES 

loop; 

AboveThreshold (loop) =0; 

BelowThreshold (loop) =0; 
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Threshold (loop) =0; 
MaxRxy (loop) =0; 
RxyThreelag (loop) =0; 


EEEKEEEEEEES Set Parameters TEEESESEEEEEEEEEEEEEEEES 

N=100; $100 noise points, therefore 200 correlation points 
a=10; %+- 10 % on either side of 0 time lag 

Pfa=PFA; $Reestablish because Pfa destroyed each loop 
PfaDefined=Pfa; 


SEES Make N points Noise $$$$55555155559559515455 
snrnoiseloop $call signal plus noise generation at a SNR 


%$%%%% Correlation of zero mean noise only for Pfa calcul. $$ 
MeanRayl=mean (RayleighNoise); $noise created in snrnoise.m 
MeanRay2=mean (RayleighNoise2); 
RayleighNoise=RayleighNoise-MeanRayl; 
RayleighNoise2=RayleighNoise2-MeanRay2; 


Rxy=xcorr (RayleighNoise, RayleighNoise2); %xcorr signal+noise 
later 


EEEEEEEEKEESEESEEEESESE Compute t~ BEEEEEEEEEELEGEEEEEEEES 
te 
for i=(N-a) :1: (N+a) 
t=tmRxy (1) ; 
end 
t=t / (28331, 


EEEEEEEEEEESEES Compute variance t EEEEEEEEEEEEEES 
Vart=0; 
for i-(N-a):1: (N*a) 
Vart=Vart+ (Rxy (1) -t) .%2; 
end 
Vart=Vart/ (2*a+l); 


EEEE%EEEEESES$SEEE$EEE Define PLA EEEEEEEEEEEEEEEEESES 

Pfa=Pfa*2; $multiply by 2 

Pfa=1-Pfa; subtract one 
Threshold (loop) =inverf (Pfa) ; 

Threshold (loop) =Threshold (loop) *sqrt (2) *sqrt (Vart); %mult by 
the sqrt (2) 


%%%%% Xcorr of signal + zero mean noise %%%%%%3%33%% 
Rxy=xcorr (DelayedPulseTrain,DelayedPulseTrain2); 
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Errr% Compare +a to -a to Threshold 2%32%%3%33%3%3%%%3%% 
for 1=(N-a):1: (N+a) 
if Rxy(1)>=Threshold(loop) 
AboveThreshold (loop) =AboveThreshold (loop) +1; 
else 
BelowThreshold (loop) =BelowThreshold (loop) +1; 
end 


if Rxy (1) >MaxRxy (loop) 
MaxRxy (loop) =Rxy (1); 
SaveThreshold (loop) =Threshold (loop); 
%SaveVart (loop) =sqrt (Vart); 
end 
end 


RxyThreelag (loop) =Rxy (97); 
end $end # Iterations loop 


$3%3%%%% Count correct threshold crossings %%%%%%%% 
CorrectTDOA=0; 
for count=1:1:Iterations 

if RxyThreelag (count) ==MaxRxy (count) 

Correct TDOA=CorrectTDOA+1; 

end 
end 
PercentCorrectTDOA= (Correct TDOA*100) /Iterations; 


diary flag4 

PercentCorrectTDOA 

BelowThreshold 

AboveThreshold 

PfaDefined 

SNR 

SaveThreshold compare threshold to zeroth 
lag below 

RxyThreelag Channel 2 lags Channel 1 by 
three pulses. 

MaxRxy 

diary off 


end $end $ test points loop 


end +end Pfa loop 
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$555555955959595459595959595959595959595955959555959595959595959595955959595959595959595959595959595595953 


% Program Name : snrnoiseloop.m 19 JULY 
EEEEEEEEEEEEE EEE EEE EEEEEEEEEEEEEEEE EEE EEE EEEEEEEEEEEEEEEEEES 
dly=0; $Channel 1 delay 

wi2=5; $Channel pulse width 


pr2=10; %Channel 1 pulse period 
nu2-10;  $Channel 1 number of pulses 
NumberPointsDel= (pr2*nu2)+dly; %calc. nmbr pts sig + delay 


dlyCh223;  $Channel 2 delay 

wiCh2=5; $Channel 2 pulse width 

prCh2=10;  $Channel 2 pulse period 

nuCh2210;  $Channel 2 number of pulses 
NumberPointsDel2=(prCh2*nuCh2)+dlyCh2; $nmbr pts sig + delay 


ZeroPadPoints=NumberPointsDel2+NumberPointsDel-1;2%N3=N2+N1-1 
for m=2:1:19; %2^(19) = 524288 point FFT max 
if 2^m >= ZeroPadPoints, 
ZeroPadPoints=2*m; 
break; $when find a 2^m power, exit loop 
end 
end 


%%%% Create Channel 1 Delayed Pulse Train %%%%%¥%%¥%%%% 
DelayedPulseTrain-zeros(1:ZeroPadPoints);$0 pad past pr2*nu2 
for PulseCounter=(dly+t1) :pr2: (pr2*nu2+(dly+t1)) ;$pr2*nu2=Nmbr 
for CutUpPulse=0:pr2; %...pts in Del. Pulse Train 
if CutUpPulse<=wi2 
DelayedPulseTrain (PulseCounter+CutUpPulse)=1.0; 
end 
end 
end 


$%%%% Create Channel 2 Delayed Pulse Train %%%%%%%%%% 
DelayedPulseTrain2-2zeros(1:ZeroPadPoints);$0 padpast pr2*nu2 
for PulseCounter= (dlyCh2+1) :prCh2: (prCh2*nuch2+ (dlyCh2+1)); 
for CutUpPulse=0:prCh2; 
if CutUpPulse<=wiCh2 
DelayedPulseTrain2 (PulseCounter+CutUpPulse) =1.0; 
end 
end 
end 


EFEEEEEEEEEEEEEEEEEEEEEE EEE EEE EE EE EEE EE EE ESEEEEEEEEEEEEEEEES 
% Signal to Noise Ratio Calculation 
% Define Rayleigh noise power to equal 2. 
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% Rician pdf created for sinusoid + Gaussian noise. 
EE EEEEEEEREEEEEEEEE EEE EE EE EEE EEE EERE EE EERE EEEEEEKEEEEKEEEEKEES 
echo off 


NoisePwr=2; $rayleigh noise pwr =2 


EEEFEFEEEEEEE EEE EEE EE EE EEE EEE EE EE EEE EE EEE EE EEEEEEEEEEEEEESES 

% Sum Channel 1 & 2 signals plus their initial delay. 

ESEEEEEEEEESE EEE EEE E EERE EEE E EEE E EEE E EEE EEE E EERE EEE EEEEEEEEES 

SigEnergy=0; $initialize Channel 1 energy to 0 

for index-1:NumberPointsDel;$sum over Ch. 1 delayed signal 
SigEnergy=SigEnergy+ ( (DelayedPulseTrain (index) ) .%2); 

end 


SigEnergy2=0; initialize Channel 2 energy to 0 
for index=1:NumberPointsDel2;%sum over Ch. 2 delayed signal 

SigEnergy2=SigEnergy2+ ( (DelayedPulseTrain2Z2 (index) ) .*2); 
end 


EKEKEEEESESEEEEEKEE EE EEE EE EE SEE EE EES EE EEE ES EE EEE EE EEE EEERES 
%* Calc. Ch. 1 signal amplitude to meet required input SNR. 
SEEEESESSESEESES EEL EE EEE EEE EEES ES ES EES EEEEEESEEEESEEEEEE EES 
SNR10=SNR/10; 

SNRLinearz10^SNR10; $SNR in linear units 
SigAmplitude=sqrt (NoisePwr*SNRLinear*NumberPointsDel/SigEner 


gy); 


EESESEESESEESEESES EE EES EE EE EES ESE EE EE EEE ES ES EEE EEESEEEESESES 
$ Calc. Ch.2 signal amplitude to meet required input SNR. 
EEEEEEEEEEEEEEES EEE ES EE EE SEE ES EE EEE EE EEE EEE GEE EEEEEEFEEEEES 
SNR10=SNR/10; 

SNRLinear=10”SNR10; $SNR in linear units 
SigAmplitude2=sqrt (NoisePwr*SNRLinear*NumberPointsDel2/SigEn 
ergy2); 


$ 35555 9,9 9 33S S€**X*X3X553559595$5$35$3535959559959593533959595959595935*99* 
% Create Rayleigh Noise matrix for Channel 1. 
EEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEEE EEE EE EEEEEEEEEEEEEEEEEEEES 
for index=1:NumberPointsDel, 

GaussNoisel (index) =rand; %N (0,1) 

GaussNoise2 (index) =rand; SN (0,1) 
end 
RayleighNoise=sqrt ( (GaussNoisel) .^2+ (GaussNoise2).^2); 


$1$51595955*54545959595995959595959595955959595955959595959559595 55 9 9 9 5999959599 99995995959599559 
$ Create Rayleigh Noise matrix for Channel 2. 
EEEEEEEEE EEL EEEEEEEEEE EEE EE EEE EE EEE EE EEEEEEEEEEEEEEEEEEEEES 
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for index=1:NumberPointsDel2, 

GaussNoise3 (index) =rand; $N (0,1) 

GaussNoise4 (index) =rand; %$N (0,1) 
end 
RayleighNoise2=sqrt ( (GaussNoise3) .*2+ (GaussNoise4) .*2); 


EEEEEEEEEEEEEE SEES EL EE EE EEE EE EE EEE EE EEEEEEE EE EE EEEEEEEELESS 
% Make Rician amplitude of Channel 1 delayed pulse train 

% and remove the dc component of entire signal. 
EEEEEEEEEEEEEE EE EE EES EEES EE EEE EEE EE EE EEE EEE EEE EE EE EEEEEEESS 
SE=0; 

for index-1:NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain (index) == *1f one, make Rician 


DelayedPulseTrain (index) =DelayedPulseTrain (index) .*rician(Si 
gAmplitude); 
SE- (DelayedPulseTrain(index)).^2-*SE; 
end $end Rician modification loop 
end 
SP=SE/NumberPointsDel; 


%$%%%%%% Add noise to Channel 1 Delayed Pulse Train %%%%%%% 
for index=1:NumberPointsDel, %1 --> delay + delayed signal 
if DelayedPulseTrain (index) ==0 %*no Rayleigh noise to signal 


DelayedPulseTrain (index) =DelayedPulseTrain (index) +RayleighNo 
ise (index) ; 

end 

end 


%$%%%% Remove dc component of Channel 1 Pulse Train %%%% 
Chldczmean (DelayedPulseTrain(1:NumberPointsDel)); 
DelayedPulseTrain(1:NumberPointsDel)-DelayedPulseTrain(1:Num 
berPointsDel)-Chldc; 


$$5959559559545555959555595459595595959595595595995959595959595959595959559595959959959959*95959 
$ Make Rician amplitude of Channel 2 delayed pulse train 

$ and remove dc component of entire signal. 
$$5559555959595955959595459559595959595*95595959595955959595*959595 959595 95959555595*59595*5559515*55 
SE=0; 

for index=1:NumberPointsDel2, $1 --> delay + delayed signal 
if DelayedPulseTrain2 (index) == $if one, make Rician 


DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) .*rician( 
SigAmplitude2); 
SE= (DelayedPulseTrain2 (index) ) .*2+SE; 
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End send Ricsan modification loop 
end 
SP-SE/NumberPointsDel2; 


%$%%%%% Add nc.se to Channel 2 Delayed Pulse Train %%%%%%%% 
for index=1:NumberPointsDel2, %1 --> delay + delayed signal 
if DelayedPulseTrain2 (index)==0 %no Rayleigh noise to sig 


DelayedPulseTrain2 (index) =DelayedPulseTrain2 (index) +Rayleigh 
Noise2 (index) ; 

end 

end 


%$%%%% Remove dc component of Channel 2 Pulse Train %%%% 
Ch2dc=mean (DelayedPulseTrain2 (1:NumberPointsDel2)); 
DelayedPulseTrain2(1:NumberPointsDel2)=DelayedPulseTrain2 (1: 
NumberPointsDel2) -Ch2dc; 


DelayedPulseTrain=DelayedPulseTrain (1:100); 
DelayedPulseTrain2=DelayedPulseTrain2 (1:100); 


RayleighNoise-RayleighNoise(1:100); 
RayleighNoise2-RayleighNoise2 (1:100); 


87 


Time Domain Correlation 


in Noise Flowchart 


| pt a | 
Detine Pta, SN Compare Axy 0) to 


30 triais per SNA | 


Threshold 


Detine Correiation 





boundary. a 







< Ru. 0 = Any m? 
mas 


= incorrect TDOA 
M estimate 
y 
Y 


Correct TDOA estimate 





^. 


| Create shifted Rayleigh 













noise 


| Axy 0) snitted Rayiergh 


Noise 


alculate noise variance 





Cakulate Threshold 
for Pta 





Create Rician delayed 
signal for specified 
SNR (delay 2 m) 







Raxy © Rician signa! 
pus shifted Rayleigh 
noise 





Figure 25. CorrPfaLoop.m MATLAB flowchart. 
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APPENDIX C. TDOA CONSTANT THRESHOLD SIMULATION 


A. INTRODUCTION 


Four graphs of the output of a time domain based 
correlation detector using a constant threshold are given. 
Each figure has four simulations using a SNR of 1dB. The P,, 


is varied from 0.01 to 0.00001 from Figure 26 to Figure 29. 


89 


de] suit) 
01 os X 0 E a 
y " 1000=8Jd "INS HP [+ 
UZ- 
Ol- 
0 
E i p8p = proysan 701 
0c 
t2 1 Mo jr. a = E 0t 
PJOYSIIYI YW uone[31102$8019 uteulop dW) 
Be] 3111) 
A OS 0 0S- _ WR. 
i | 100 —UJd "NS PL 
E 07- 
| 01- 
0 
Ops? F Frsamqp j01 
-107 
Ot 
| | Up 


P[0YS231Y) YJIM UONEL3110955049 UNLWOP IUN 


(DAY 


(DAY 








avy su) 
000 0S 0 Us; u 
" 10'0=8Jd “UNS HPI 
0c- 
Ol- 
[i 0 8 
| = 
PE ux per AGILE HE DIO 
TTT | pi 
-107 
— d us i SS 9. ——— 0t 
PJOYSIIYy) YW uonegj3109$$012 Uteulop ault 
de| aui 
06 0 UA Ope- 
| " 100=8d UNS PI + 
07° 
Ol- 
10 


(Ary 


ddr & proysaryy 10! 


02 


30£ 


= i 


l ' l z Or 
pI[ousa1u UYJIM UONELIJIOISSOID uteulop un 


Time domain based correlation detector output 


for a SNR of 1dB and a Pr or 


Figure 26. 


90 


se Iun 


001 Us 0 05° 001- 001 
| De LLLI —— -0b- 
100 0=Ud “UNS HPI + 
0c 
Y wi 
AA d PON OEA PL - pioussuu - 
/ i pjousaa L m 
I t EA a OP 
pIousa2u) ui uone[341102580J2 UIEWOP Just) 
avy aw 
001 Us Ü DS- DT 001 
— = T 5 B aM ES Zar 


'L00'0 =8d UNS UPI+ 





— 
PLOYSIIY) YUM uon [3410555042 ureulop aui Ut 





Jel awn 


OS 0 os- 


—q[-————- —-—p 


IOWO =8V]d UNS PI + 


L Som eh ze 


: A ET | Er 
pIoysasyı yım UONEL31I0955039 wewop awn 
de au) 


u 29 O 
E "1000 =8Jd ANS YPI + 


2841 = PIOYS3ML 


L : l i = 
pjoysaay) Y)IM UONLL3II0985019 UNEWOP Iwn 


001- 


üc- 


09 


“be: 
07- 
ot- 

40 

ol 
07 
0€ 


Ov 


MAy 


MAy 


detector output 


Ton 


* 


domain based correlat 


Time 


Figure 27. 


Of 0 UO 


hor aMSNR of IdB and a P; 


9T 


dv| aut) 


001 NU 0 0s- 


E 





otolla uns pre [07 
0 
VA deed E prono 
BE Joz 
Fer A A 20 AAA en OP 
pfoysaayı uit UOHe]ILIOISSOIN UlRWOp UIT) 
de] auım 
ool o5 0 05- 


RC. 001 
1000°9=8d “UNS HPI * 


= m 01- 


301 


BLOT = plogsary 
Y Lae. 


l l : L S 
pjoysasyı yım uon[3J102885012 Ureuiop oui vi 


(Ary 


(DAY 


001 os 0 os- 


dv] aun 
001- 


— —- jt 


——-- —— + 


Ioood=bia uns apı+ 70° 


OT- 


| 401 
as WV | egz F Proussnur - 


0c 





TE SH ee L - = ¡AE 0€ 
pIousa1u) ui^ UOHn?[231102$8012 utetuop aui 


Be] awn 


€ 090  — 0€ OK. 


———T—— ¡— 


odo =e1a uns apı+ [OO 
Ot 
10 
ol 


YW YW ice br = prussn- 
0c 


0t 


1 " 1 wl ES OP 
pioysasy uU UOI ?[31102855012 Uureuiop aui 


(Ary 


MAy 


Time domain based correlation detector output 
of 020001 


for a SNR of 1dB and a Pz, 


Figure 28. 


92 


dv] aun 


A EE US NL 001 OS 
| I(o000 2o1d ^3NS api [Y E S2 

oz- 

0l- 


(DAXA 


-OI 


6 = pyoysamyy 10€ 


de) au 


0 os- 001- 
——p—— 4p —— qot 


100000 ="Jd “UNS API + 
: OZ- 


| | LE6S RE = PIOYSIMUL 02 


—— n es d -—]H -—— — AS ———— en OL Ec EET. age ee o UE E E OF 
pjoysaayı yım uone[341102$50J2 uteulop ouln piousa1q) qii uon e[34102858s019 Uleulop aul) 
sep Iwn dey ow 
NAMEN uM AME OS 100 a Le EB END BE p. 
! À i j 100000 =8]d ANS HPI+ 
| Ob- 107- 
rg JA UNS HPI + 
-JOc- or 
/ 0 $ 1 0 
^| | [DHS E 107 OT 
OP BPLCGT — ppusanjp 710€ 
| +09 D€ 


= l 
PIOYS31Y) YIUM uonpp 31102858022 ureiuop aug 


e zu ; i ps 
poys YlM VONPJIJIOISSOJO Utputrop oum] 


(DAY 


(DAY 


detector output 


KON 


domain based correlat 


mereamoiie of IdB and a P, of 0.00001. 


Time 


Figure 29. 


95 


APPENDIX D. FFT OUTPUT NOISE TO SIGNAL RATIO 


A. INTRODUCTION 

The noise to signal ratio for thirteen passes through a 
FFT butterfly are calculated below. The calculation assumes 
the word length of the arithmetic logic unit (ALU) of the 
computer is 16 bits long. The length of the input data points 
are 12 bits long. The data points are right justified when 
placed into the ALU. For a 16 bit ALU and a 12 bit data 
input, scaling occurs at the seventh pass and every other pass 
after that. 
Pass 1l 

The noise power at the output of the butterfly after the 


first pass N, (1) is a function of 


Ni(1) = 2N?(0) + 202(0) + S202(0) + 0%(0) 


where S* = signal power input to butterfly 
o> = truncation noise power (74) 


o% = scaling noise power 


oĉ = sine/cosine table noise power 


N (0) = noise power input to FFT 


During pass one there is no scaling, no truncation, and the 
trig table is not used. These error terms are zero. Equation 


74 reduces to 
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Na(1) = 20, 
^ mes 
where N, = a =) 


noise input to butterfly 


Pass 2 


The noise power at the output of the butterfly after the 


second pass N,*(2) is a function of 


Ni(2) -2NSU - 202 (07) 255?*020" 020) . (76) 
During pass two scaling and truncation are not performed, and 


the trig table is not used. These noise terms are zero. 


Substituting Equation 75 into 76 


Na (2) 


ZNE) bis 


DENE 


Pass 3 


The noise power at the output of the butterfly after the 


third pass N,*(3) is a function of 


[89 2N:(2) *205(2) + 2250712) + 0712) . (78) 
Scaling noise is zero because it is not performed during pass 


three. Substituting Equation 77 into 78 


35 


Ni(3) - 2(2?N2) + zog ea, 


(79) 


SNE ada sto NEM 


Pass 4 


The noise power at the output of the butterfly after the 


fourth pass N,/(4) is a function of 


Nå (4) = 2N (3) + 202(3) - 2? s?o (oo tee (80) 
Scaling noise is zero because it is not performed during pass 


four. Substituting Equationz79zınressy 


Nx (4) - 2*NZ *"295?g. ES TN (81) 


Pass 5 


The noise power at the output of the butterfly after the 


fifth pass N,/(5) is a function of 


N2(5) - 2Nf(4) + 202(4) + 24520%(4) + o*(4) . (82) 
Scaling noise is zero because it is not performed during pass 


five. Substituting Equation 81 into 82 


NA(S) 9 23N. v (2999 2: y Seq OP (83) 


Pass 6 
The noise power at the output of the butterfly after the 
sixth pass N,'(6) is a function of 
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DNE ME o. 1145.20 5) £955) . (84) 
Scaling noise is zero because it is not performed during pass 


Si ubstitutzrng Equation 83»into 84 


Na 


= 2 2 2 
le) 2. N 227520, t+ 1507... (85) 
Pass 7 

The noise power at the output of the butterfly after the 


seventh pass N,’(7) is a function of 


N3(7) = 2N2 (6) + 202(6) + 26520%(6) + 0%(6) . (86) 
Scaling is performed for the first time during pass seven. 


Substituting Equation 85 into 86 


Be Dun (285: 26,9520%.+ 3107 + 205 . (87) 


Pass 8 
The noise power at the output of the butterfly after the 


eighth pass N,/(8) is a function of 
NZ(8) = 2NZ{(7) + 202(7) + 27S202(7) * 0207) . (88) 
Scaling is not performed during this pass. Only the scaling 


noise from the previous pass is added. Substituting Equation 


87 into 88 


NÈ (8) = 29N + (2° + 2°) S?o% + 6307 + 405 . (89) 
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Pass 9 
The noise power at the output of the butterfly after the 


ninth pass N,/(9) is a function of 


(8) * 202(8) + 28520%(8) + 05(8) . (90) 


Substituting Equation 89 intorg 


NA (9) = 2903 + (21° + 29 + 2%) 920% + 12707, + 10g, NE 


Pass 10 
The noise power at the output of the butterfly after the 


tenth pass N,’(10) is a function of 


N2(10) = 202 (9) + 202(9) + 2%5202(9) + 02(9) . (92) 
Scaling is not performed this pass. Only the scaling noise 


from the previous pass is added. Substituting Equation 91 into 


92 


Ng(10) = 210, + 212520, + 2590 00 (93) 
Equation 93 describes the noise power at the output of a 
butterfly after the tenth data pass. The noise to signal ratio 
after the tenth pass through the butterfly is expressed as the 


ratio 
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N? Bene 27550, 1200: 
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S? 210g? 


Dividing through by the denominator (signal power) each 


independent and uncorrelated noise term can be identified 


2 2 2 2 
Na (10) Mc + 40% + Beer + 1 Is 
Ss? S? 4 s2 IES 
Nx | | 
where E (10) = output Noise to Signal ratio 
NS 
— = input Noise to Signal ratio (95) 
Ss 
40% = sine/cosine table noise 
i or , : : 
s -truncation Noise to Signal ratio 
2 
1 Os _ 5 1 ; 
3192 -scaling Noise to Signal ratio 


The noise to signal ratio calculation for 13 passes (2^ 


= 8192 points) through the butterfly algorithm continues 


Pass 11 


The noise power at the output of the butterfly after the 


11% pass N (11) is a function of 


N: (11) = 2N2 (10) + 202(10) + 21%520%(10) + 0(10) . (96) 


Substituting Equation 93 into 96 


DE ati (2)? 219 529.2 5119» 5* 420. . (97) 
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Pass 12 
The noise power at the output of the butterfly after the 
12% pass N, (12) is azfunet vongor 
Ni (12) = 2N/(11) + 20%3(11) + 2t!s204(11) + 02(11) . (98) 
Scaling is not performed this pass. Only scaling noise from 


the previous pass is added. Substituting Equation 97 into 98 


N, (12) = 212N, + (220 + 212) 520, -= 10200 war ve (99) 


Pass 13 
The noise power at the output of the butterfly after the 


13°" pass N, (13) is a function of 


(100) 
Ni (13) = 2Ni(12) + 20%(12) + 212529%(12) + 02(12) 


Substituting Equation 99 into 100 


(101) 
Ni (13) = 25 + (215 « 219 + 212) sig 20410 100 H 


Equation 101 describes the noise power at the output of a 
butterfly after the 13'^ data pass. The noise to signal ratio 


after the 13'^ pass is expressed as the ratio 


Na _ 2NZ + (235 + 29 + 212) 570), + 204707 + 1700: 102) 
Nep A 21752 l 


Dividing through by the denominator, each independent and 


uncorrelated noise term can be identified 


100 


1 Os 
48.2 S? 


z 2 

1 Or 1 Os 
O tem m 
er 52 48.2 52 


IN 


| Rt 


output Noise to Signal ratio 


- input Noise to Signal ratio 


- sine/cosine table noise 


-truncation Noise to Signal ratio 


=scaling Noise to Signal ratio 


DOT 


(103) 
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